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Abstract 


We have carried out ab initio total energy calculations using the density-functional theory 
to study the native defects in GaAs material. First principle pseudopotentials for Ga and 
As were generated according to the Bachlet-Hamann-Schliiter scheme using the computer 
code fhi98pp. These pseudopotentials were fully separable and norm conserving in 
nature. In order to ensure the accuracy of the results hard pseudopotentials were 
generated. The total energy calculation was performed with another code flii96md. 
Lattice constants of Ga, As and GaAs were optimized on the basis of total energy 
minimization. All the calculations were done at 20 Ry cutoff and only the electron 
minimization scheme was employed using the second order damped Joannopoulos 
algorithm. The supercells for the calculations were constructed from four cells of Ga and 
As each and eight cells of GaAs. For exchange and correlation local density 
approximation (LD A) was employed. Heat of formation of GaAs was calculated for these 
optimized structures by subtracting the bulk energies of Ga and As from the bulk energy 
of GaAs. Though the individual energies of the materials concerned are close to those 
published in the literature, our results show that the heat of formation is 0.36 eV where as 
the experimental result is 0.74 eV. Various point defects were introduced in a GaAs 
crystal by physically removing the atoms. The formation energies for defects at different 
charge states have also been examined. 
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Chapter 1 

Introduction 

1.1 Background 


The calculation of total energy using the density functional theory (DFT) [1.1] is 
the most commonly used tool in condensed matter physics and materials science. This 
method has been successfully applied to molecules, bulk solids and surfaces and the 
results are in very good agreement with the experiments. The basic essence of this theory 
is that, given an external potential, we are in a position to determine the electronic 
density, energy and any ground state property. Electronic structure theory, i.e., the 
theoretical prediction of materials properties from fundamental quantum mechanics 
without recourse to empirical parameters, has advanced significantly in recent years, 
serving basic science as well as applied research in numerous industnes. The 
applicability of current methods and their spin-offs covers an extraordinarily wide range 
of fields and interdisciplinary research topics that embrace Condensed Matter Physics, 
Materials Science, Computer Science, Chemistry, and Biology. 

Defects in semiconductor are particularly important in the context of optical and 
electrical conductivity. As an example, the semi-insulating GaAs has applications m 
hetrostracture bipolar transistors and optical devices. 

In ‘’ab-inito" or ‘first-principle’ calculations an exact description of the system is 
used i.e. no empirical parameters are used. In atomic ab-initio calculations only some 
basic data such as atomic numbers and masses are used as inputs. The ground state of the 
system is then obtained from DFT using the Bom-Oppenheimer approximation [1.2]. The 
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corresponding equation describing the system is known as the Kohn-Sham equation [1.3]. 
The ground state energy is then obtained by minimizing the energy functional. The DFT 
formulation of the Kohn-Sham equation is an exact one. However, the exact form of the 
exchange-correlation energy as appearing in this equation is not known. Certain physical 
approximations are made to approximate these energies. A widely used approximation is 
the local density approximation (LDA). Though LDA underestimates the band gap 
energies of the semiconductors, it accurately predicts the ground state geometnes, charge 
densities and most importantly, gives the physical trends correctly. Beyond LDA, an 
alternative method, generalized gradient approximation (GGA) [1.4] was presented, 
which gives accurate results for the semiconductor systems. But these GGA calculations 
are computationally more intensive. 

1.2 Literature Review 

k 

A good description of the fhi96md code, which we use here, was given by 

j 

Bockstedte et. a/.[1.5]. The formulation and the generation of the corresponding 
pseudopotentials xi^,mgfhi98pp code was described by Fuchs et. al. [1.6]. 

Chadi and others [1.7] had calculated the total energy and the equilibrium lattice 
constant for Ga, As and GaAs in a different context (surface reconstruction of the GaAs 
(100) surface). They had also calculated the heat of formation as a check for their 
optimized structure. The chemical potential, jj. is defined as the derivative of the Gibbs 
free energy, G = £ + FV - TS , with respect to the number of particles. Here the term 
PV is negligible and at temperature T = 0, G = E . The heat of formation of GaAs is 
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given by, ^ = MGaAs(buik) - Moaibuit) -M^buik)- Since LDA gives accurate results when 

the energy differences are concerned, they had calculated the total energy of the 
concerned materials independently and used the difference of the energies to get the heat 
of formation. 

They had observed that though this kind of calculations are mainly done at a 
cutoff of 12-18 Ry. But since, the surface contains large number of atoms, which for 
equivalent number of plane waves requires a smaller cutoff energy, surface related 
calculations could be done at 6 Ry. Instead of doing the surface calculation at 6 Ry cutoff 
with hard pseudopotentials, Chadi and others had carried out the calculations with softer 
pseudopotentials, which converges at 6 Ry cutoff. They had modified the earlier Kerker 
approach [1.8] by choosing a larger cutoff radius and thus producing a softer 
pseudopotential. The pseudopotentials still are norm conserving. In order to ensure the 
accuracy of this softer pseudo-potential, various bulk properties were calculated. The 
GaAs stracture and related properties were calculated first. The calculation was done 
with 8 special k-points. The results show that when a hard pseudopotential predicts a 
lattice constant of 5.42 A and 5.55 A at 6 and 12 Ry cutoff, respectively, the same is 5.48 
A and 5.52 A for the corresponding softer one. The phonon frequencies calculated at a 
higher cutoff of 18 Ry had a difference of only 2.5%. Similarly the heat of formation for 
these two cases were compared and the results are fairly in good agreement with 
experiment. This will be discussed in a greater detail in chapter 4. 

The As stracture related calculations were done according to the method 
prescribed by Nielsen and others [1.9]. According to Nielsen, the ground state of As is 
A7 structure. The structure may also be described according to simple cubic stracture, but 


3 



this structure is not a stable one. Chadi had calculated the As total energy using the 
experimental structure with 44 inequivalent k points. 

The structure of Ga in its ground state is orthorhombic in nature. The total energy 
for this structure is calculated using 75 inequivalent k points at a cutoff of 6 Ry. 


1.3 Objectives of the present study 


The objective of the present work is to estimate the various point defect formation 
energies in GaAs. For this purpose, a norm conserving fully separable pseudopotentials 
for Ga and As are generated using the computer code fhi98pp. The pseudopotentials and 
the pseudo wavefunctions are input to another code fhi96md for calculating the total 
energies of the Ga, As and GaAs. The structures are to be optimized on the basis of total 
energy calculation. The accuracy of the structures is verified by calculating the heat of 
formation. Then the various defects of GaAs are studied. 
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Chapter 2 

Electronic structure calculation 

2.1 The problem of the structure of matter 

The microscopic description of the physical and chemical properties of a material 
system is a complex problem. In general, these systems involve a collection of atoms 
interacting with each other. These systems may be isolated, extended or a combination of 
both. The Hamiltonian of such a system may be written as, 


H= -y—vj -y- . ... _ 


P P 7 7 

-v? 


^2 N N 


+ 


N N % P N 

7 ~\ 


2 ; 


2 - r 


/~1 /*=! 


( 2 . 1 ) 


where R={Rj},I =1,2, P, is a set of P nuclear coordinates, F = {/;},/ = 1,2, N , is a 

set of N electronic coordinates and e is the electronic charge. Zj and M, are the charge 
and mass of the / th nucleon while m is the mass of the electron. In principle, all the 
properties of the system ean be derived if we can solve the following Schrodinger 
equation, 

H'T(f,R)=E'F(f,R) (2.2) 


where E being the eigenvalue corresponding to the eigenfunction 'F(f,R). In practice, 
this problem is impossible to solve in the full quantum mechanical framework. And there 
are very few possible cases where the above equation can be solved analytically or 
numerically. Common is solution in the ground state of such a problem under certain 
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approximations. The usual choice is to use some realistic approximation. 


2.1.1 Adiabatic approximation (Born-Oppenheimer approximation) 

The first observation is that all the nuclei, being heavier than the electrons, move 
much slower than the associated electrons. As a result, at any instant, the nucleus may be 
considered to be at rest while the electrons are rotating around them. In other words, as 
the nuclei follow their dynamics, the electrons associated with them adjust their 
respective wave functions according to the nuclear wave function. If we substitute, 

(2.3) 


where ©„(.R,t) is the nuclear and <^„(r,R) is the electronic eigenfunction, the electronic 

and nuclear part separates out. This approximation is valid for m/M or . The 
corresponding electronic Hamiltonian can be written as. 


N N N -I P N y 

,=[ 2m 2 ,=i j^, r, — /=i ,=i R^ — r, 


(2.4) 


and the Schrddinger equation reads as, 

He(I>(?,R)=Ee<l>(f,R) (2.5) 

where Eg is the energy eigenvalue associated with the electron. The above 

approximation is known as the adiabatic approximation or the Bom-Oppenheimer 
approximation [1.2]. 


2.1.2 Quantum many body theory 

In order to solve Eq. 2.5, the first approximation method was proposed by Hartree 
[2.1]. According to this theory, the many body wave function can be written as a product 
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of one-electron wave functions i.e. 


In this theory, the atom-electron interaction is replaced by a effective potential. 


V +V^j.{r,R)\(pXr,R) = £,^,{r,R) 


where 


= (r,R) + f , 

^ \r - r \ 


p{r)=^Xr,Ry 


is the electronic density associated with the j th electron. The total energy of the system 
is given by, 




( 2 . 10 ) 


The above method is known as self-consistent Hartree (SCH) method. The main problem 
with the Hartree theory is that it treats the electrons as distinguishable particles. The 
indistinguisability of the electrons was taken into account in Hartree-Fock (HF) theory 
[2.2, 2.3]. Here, instead of writing the many electron wave function as a simple product, 
the total wave function is written in the form of Slater determinant. 


^r,R)=SD{^XrM- 


(piiR^rX (PiiR^r^) ••• 
]_<Pzi^rX <pXKr2) ••• 

i i 


( 2 . 11 ) 
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However, the Hartee-Fock theory totally ignores the Coulomb correlation, which arises 
due to the repulsion between the electrons. The most straightforward way to introduce the 
electronic repulsion is to write the many body wave function as a linear combination of 
the Slater determinants ( i.e. Hartree-Fock wave functions). This method is known as 
configuration-interaction method (Cl). 

2.1.3 Density Functional Theory 

Parallel to the development of this Hartree and related theories, Thomas and 
Fermi [2.4, 2.5], at the same time as of Hartree, proposed that the full electronic density 
is the fundamental variable, instead of the wave function. They derived a different 
differential equation for the density without resorting to the one-electron orbitals. A firm 
and exact theoretical formulation of such description was provided by Hohenberg and 
Kohn [1.1] in 1964 and is know as the Density Functional Theory (DFT). In 1965, Kohn 
and Sham [1.3] proposed the idea of replacing the kinetic energy term in Eq. 2.5 by an 
equivalent of a non-interacting system. Kohn-Sham defined the density p{r) as, 

where is the occupation number of the orbital (Z?, ^ (F) . Such a procedure brings us to 
self-consistent Kohn-Sham equations: 

(^.s(r)=s,_,^,JF) (2.13) 

where the effective potential y[p; r] is written as. 
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(2.14) 


np-.n=V„(r)+ \,^d'r'+v,^[p] 

^\r — r\ 

=-^+yM,n+y,c\p(ri\ 

(p^ ,{r) are known as Kohn-Sham orbitals and they may or may not have any physical 
significance. is the functional derivative of the exchange and correlation energy, 
Exc ■ The total energy of the system is then given by, 

jyMP(r)d’r (2.15) 

The problem with Eq. 2.15 is that the exact form of E^c (or V^c) is not known. Some 

physical approximations are made to calculate these quantities. These approximations are 
known as Local Density Approximation (LDA) and are very successful in electronic 
structure calculation. Using LDA, V^c is written as, 

dp{r) dp 

is known as the exchange and correlation energy density and it has been 

calculated by using quantum Monte Carlo method by Ceperley and Alder [2.6] and 
parameterized by Perdew and Zunger [2.7]. These are considered to be most accurate to 
date and are being used in our calculation. 

2.1.4 Dynamics of ions 

Once the ground state of the electrons are calculated as described in the above 
section, the atomic equation of motion is integrated with standard molecular dynamics 
(MD) techniques. The choice of the scheme depends on the type of the application. 
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Among various MD methods, the most popular ones are Verlet [2.8] and predictor- 
corrector [2.9] method. 


2.2 The concept of pseudopotential 


The electron-nuclear interaction in an atom is defined by the Coulomb interaction. 
The electrons in an atom can be classified into two classes: one of them moves freely 
through the nucleus and are responsible for the chemical properties of the matter i.e. the 
valance electrons, while the others are tightly bound to the nucleus and are called the core 
electrons. In a calculation where both of them are taken into account is called all-electron 
calculation. Since the core electrons do not take part in the chemical bonding i.e. does not 
guide the chemical properties of the atom it is possible to integrate out the corresponding 
effects of the core electrons as a screening effect between the valance electron and the 
ionic cores i.e. the nucleus and the core electrons. This type calculation is called the 
pseudopotential calculation. 

The valance states, due to the orthogonalization of the core states of same 
symmetty, show remarkable oscillatory behavior with a number of nodes equal to the 
(n-.^-l), where ‘n’ is the principal quantum number and ’ is the angular quantum 
number. On the other hand, though the nodeless wavefunctions are non-oscillatory, due 
to the lack of orthogonality, they can approach the nucleus easily and as a result they 
create strongly bound states near the nucleus, which are sharply peaked close to the 
nucleus. These features are totally undesired from the point of view of the computational 
calculations, as a large number of plane waves (basis functions) are required to describe 
these noded wavefunctions, which in turn translates into vast amount of computer time 
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(the matrix to be diagonalized is too large). 


Thus the pseudopotential theory is constructed in two parts: 

• Core electrons are removed from the calculation and the interaction of the valance 
electrons with the nucleus plus the core electrons (including orthogonalization) is 
replaced by an effective, screened potential. 

• The full ion-electron interaction, which includes the orthogonality of the valance 
wavefunction with the core states, is replaced by a weaker pseudopotential that 
acts on pseudo wavefunctions rather than the true wavefunction. 

It is quite clear from the above discussion that the constraction of pseudopotential 
can be done arbitrarily. As a matter of fact, the construction of the empirical pseudo- 
potentials have been obtained in past by fitting with experimental values. But these 
pseudopotentials lacked a very fundamental property, namely that of transferability. The 
first transferable pseudopotential was constructed by Hamann, Schiiter and Chiang 
[2.10]. The modem pseudopotentials are constracted according to the recipe prescribed 
by them. These are as follows: 

1 . The pseudo wavefunctions should be node less and it should be identical to the all 
electron wavefunction outside a suitably chosen cut-off radius, re'. 

<Pps{'>') = 9psir) forr{r^ 

and 

<Ppsir) = <Paeir) forr>r^ 

2. The first and second derivatives of the pseudo wavefunction are continuous at ic. 

3 . The eigenvalues of the pseudo wavefunction should be same as that of the all- 
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electron wavefunction. 


4. The norm of the true and the pseudo wavefunctions inside the core radius should 
be same. 

5. Other conditions for the smoothness of the pseudo wavefunctions. 

Several schemes were proposed to generate first-principle pseudopotentials. Among them 
the most popular ones are Bachelet, Hamann and Schiiter [2.11], Troullier-Martins [2.12] 
and Vanderbilt [2.13]. 

The Bachelet, Hamann and Schiiter type of pseudopotential is the most 
common type of pseudopotential and are being used in the present calculation. • 

2.2.1 General background and formulation 

In generating a good pseudopotential for a material the following considerations 
are taken into account, (i) Transferability: the transferability of a pseudopotential is its 
ability to describe the valance electron configuration when placed in different 
environments. This is a major difference with the earlier empirical pseudopotentials. 
These empirical pscudopotcntials were designed for some specific environments and 
were inevitably very less or almost not transferable to different environments, (ii) 
Efficiency: one other important issue is the efficiency of the pseudopotential i.e. how it 
reduces the computational load. That means the lesser the number of basis functions it 
takes to compute the wavefunction and the charge densities, the more efficient it is. 

The transferability of a pseudopotential depends on (a) correct scattering 
properties, (b) the choice of the cutoff radius r^, (c) the adequate account of the 
exchange-correlation between the core and the valance states, (d) the Kleinman-Bylander 
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form [2.15] of non-local pseudopotential, (e) the frozen core approximation and (f) the 
treatment of higher angular moment components. The choice of cutoff radius is a crucial 
factor. The increase in the value of produces a softer pseudopotential, which converges 
more rapidly. As a result this becomes less transferable as the pseudo wavefunction, at 
the bonding radius, becomes less accurate. On the other hand should be greater than 

the outermost node of the corresponding all electron calculation. Choosing close to 
that node produces bad pseudopotential. 

2.2.2 Atomic all electron calculation 

The present pseudopotentials are constructed from atomic all electron calculation. 
The all electron wavefunctions are obtained by self-consistently solving the Kohn-Sham 
Eq. 1.3 ( here only the radial part is considered as the solution of the angular part is 
trivial), 

-l:4T+^^^+np-Ar^{r) = s„,r<^{r) (2.17) 

2 dr 2r 

where V{p',r'\ is the one electron potential. 

; r] = — + V„ [p-,r] + ipir)) (2.18) 

r 

Z 

Here is the well known coulomb interaction, /?(r) is the sum of the electronic charge 

r 

densities of the occupied states, T^[/7;r] is the Hartree potential and the 

local density approximation for the exchange-correlation potential. The total energy of 
the system is given by, 

f r+EM- (2.19) 

, 2 \r-r\ ■’ 
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= Up\+E„{p]^E,c[p\+ \--p{r)d^ 

J y. 


( 2 . 20 ) 


where Tq[p] is the kinetic energy of the non-interacting electrons, 

Exc\-P'^ is the exchange and correlation energy, 

is the electrostatic or Hartree energy. 

2 ^ 

The ground state is obtained by minimizing the total energy with respect to the density 


functional p{r) i.e. 


subject to, ^p{r)d^r = N . 


This equation gives. 


d{E[p]-N} 

dp(r) 


( 2 . 21 ) 


■^^+VAp;r]-s, <pXr)=0 


( 2 . 22 ) 


where V^[p;r]= [ ■ .— + is the effective potential. 

•’|r-r'| dp(r) 

p(r)=^n, |<z^, ()^)f , is known as the Kohn-Sham orbital which may or may not 

I 

have any physical significance. In order to solve the above secular Eq. 2.22, a spherically 
symmetric electronic charge distribution is assumed. As a result the Kohn-Sham orbitals 
separates out as. 


jR., {s,;r) 
r 


(2.23) 


where {0,(p) is the spherical harmonics. The solution of the angular part relating to 
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{9, (p) is straightforward and will be dropped in the later discussion. The radial 
equation, with the relativistic correction, relating ^ becomes. 


1 


1 dV{r) d 4£ + l) 


r 1-- 


2M{r) 1 dr 2M(r)c dr dr 


\+V{r)-£, 


(£-,;r)=0 (2.24) 


is -V(r)) 1 1 

where the relativistic electron mass is given by, — r-24 with — = 

2c c 137.036 

is the fine structure constant. The exchange-correlation potential is derived from the 
exchange-correlation energy density using 

ls!f'‘lfi]p(r)dr (2.25) 


The exchange-correlation energy density can be splitted into two parts. The exchange 
density term has been calculated by Dirac [2.14] and is given by. 


^LDA 


[/>]=-T 






,1/3 


4 ^ 4 ;^^ 


^ 1/3 


0.458 


(OM.) 


(2.26) 


r^ is the average radius containing one electron, p ' = (4.;rr/ /3) . For the correlation 

part, we have used the Ceperley and Alder [2.6] approximation parameterized by Pedrew 
and lunger [2.7]. 


,LDA 




■0.1432 


1 + 1 . 05397 ^ + 0.3334r, 
- 0.0480 +0.03 11 In r, 

+ 0.0020r Inr 


forr^ >1 

•0.0116r, 

forr^<l 


(2.27) 


2.2.3 Screened norm-conserving pseudopotential 

Once the all electron potential and wavefunctions are obtained by solving Eq. 3.7, 
the pseudopotential is constructed from it. The radial pseudo wavefiinction is then 
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2.2.4 Construction of Bachlet-Hamann-Schliiter pseudopotential 

The Bachlet-Hamann-Schliiter norm conserving pseudopotentials are constructed in 
five steps : 

(i) Kohn-Sham equation, as corrected by Dirac for relativistic effects, is solved for a 
given atomic configuration v , 



2r' ^ 




(2.30) 


/V 

(ii) Then the first step pseudopotential FJis obtained by cutting of the singularity at 
origin from the all electron potential. 

is a smooth cut-off function, which approaches 0 as 

rlr^ ->oo, cutting off around r « and approaches 1 as r->0. The optimum 
value of the parameter ^is taken as 3.5. The constant cj is adjusted to give a node 
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less wavefiinction with the same eigenvalue as that of the all electron wave- 

function. The optimum value of r^j is determined by scaling down from the 
radius of the outer most peak of the radial wavefunction G^, r, =r^,/cc, 
where cc is in the range 1.5-2. 

(ill) Once the first steps pseudo wavefunction is obtained, the second step involves 
modification of the same for the norm conservation correction for r> r . The second 
step pseudo wavefunction is given by, 

' (2.32) 


with(^j)^ + S''j f {r ! dr = \ , being the normalization 

condition from which the parameter <5J is to be determined. 

(iv) In the third step the final screened pseudopotentials V^j (r) is obtained by inverting 
the Schrodinger equation for the node less eigenfunction co^j (r) with the eigenvalue 


'' 2<(r) 


— — 2 1 — 5lL +2s, -2V;/r) 


(2.33) 


(v) Finally the screened pseudopotentials Fjj ('") unscreened using the node less 


pseudo wavefunction (f) • 




Sp^r) 


(2.34) 
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The total ionic pseudopotential is the sum of the ionic part and the spin-orbit 
coupling term as given by, 

KT (o = El (oish 1 (2-36) 

I 

where the spin-orbit term is, 

(2.37) 

For a high precision fit, the Bachlet-Hamann-Schliiter pseudopotential is split up 
into a long-range( (r ) , £ independent) and a short-range ( AF/°" (r ) , £ 
dependent) term. The core part is given by, 

Ecrer/[«"™)''V] (2.38) 

r ^ 


where c"™ = 1 . And the ionic part of the same is given by, 

Afrir) = 't{A, + A,^,)e"''' (2.39) 

J = 1 

The parameters and A, for different elements are calculated and tabulated in 
the literature. 
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2.3 Point defects in semiconductors 


There are mainly three different types of defects present in a semiconductor. The 
most common type is the vacancy defect, which is formed when an atom is missing from 
an otherwise perfect crystal. Interstitial defects are those defects in which a host atom is 
present in the interstitial site. In addition, the antisite defects are present in AB type of 
semiconductors, where an A type of atom is present at the site of atom B and vice-versa. 

Point defects in semiconductor act as a trap to the electrons and holes. As a 
result, these defects can be in different charge states. The energy required for the 
transition from one charge state to other can be studied by means of total- energy 
calculation. Here, the Fermi level acts as the reservoir. 

In actual calculation, when a defect is inserted in an otherwise perfect crystal, by 
means of periodic boundary condition of Bloch’s theorem, the defect is translated to an 
uniform array of repeated defects. This condition necessitates the use of large supercells 
in this type of calculation. 

The defect formation energy is defined as the energy required to create a defect in 
a perfect lattice. The defect formation energy in a GaAs is given by, 

, Mca > J + Qd (^v + AJ - «(?a Ag« -^As/^As (2-40) 

where and n^ are the number of Ga and As atoms in the defect cell and // ’s are 
their respective chemical potentials. is the chemical potential of the electrons ( i.e. the 
Fermi energy) and is the total energy of a supercell containing one defect at charge 
state (2o . Using, jUg^ = , which is valid in equilibrium condition, equation 

(2.40) may be written as. 
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(2.41) 


where, 

~ -^D ""^(^Ga '^^As)MGaAs{bulk) ““(^Ga ~ ^ As^i^Gaibulk) ~ ^Asibulk)) (2-42) 

~ ^Ga ~ f^As ) — i/^Ga(bulk) ~ /^As(bu!k) ) (2-43) 

E'lj values have been calculated by Northrup and Zhang [2.16]. Here, all the total 
energies are calculated about a common reference energy level. In supercell calculations, 
no absolute reference of energy or eigenvalues is present [2.17]. But, the common 
practice is to chose the top of the valance band maximum as the point of reference. 

The defect related total energy calculations, involving the periodic boundary 
conditions, needs special attention. This is due to the fact that, an uncompensated charge 
in this type of calculation, will be repeated in every possible direction and would lead to 
misleading results. To avoid this, an uniform background charge of opposite sign is 
generally supplied in the background. 
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Chapter 3 

The fhi96md and fhi98pp code 
3.1 The fhi96md program structure 


The program structure of the computer code flii96md is given in figure 3.1. The 
different blocks are self-explanatory. The basic idea of the program is to solve the Kohn- 
Sham equation self-consistently, which is described in the left hand side. The right hand 
side deals with the dynamics of the atoms. 

In the initialization part, the input files i.e. inp.mod, start. inp and pseudopotential 
files are read and the initializing input file inp.ini is being generated. Then the structure 
factors and phase factors are calculated. The initial electron density and both local and 
non-local contribution to the energy, potential and forces are calculated thereafter. Now, 
if the atoms are allowed to move, the atoms are moved according to structural 
optimization or molecular dynamics schemes. The electronic wave functions are 
calculated self-consistently, satisfying energy (epsel) or density (epsekinc) convergence 
criterion, for this new positions of the atoms. The calculation is repeated for different k- 
points. The charge density and forces are calculated once again. If the force convergence 
criterion is not met, the atoms are moved again and procedure is repeated till the criterion 
is satisfied. The program is terminated on convergence or if the number of iterations or 
cputime is exceeded. However, if termination occurs without convergence, we change the 
convergence criteria until it is satisfied. A very good description of the fhi96md code is 
given by Bockstedte and others [1.5]. 
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Fig.3.1. Flow chart of the program flii96md ( from Bockstedte et.al. [1.5]) 
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Solving the Kohn-Sham equation: Solving the Kohn-Sham equation is the most important 
part in this t 3 q)e of computer codes. The essential idea is to minimize the energy 
functional with respect to the wave function. At first, a starting guess is made for the 
initial choice of the wave functions. Then the energy minimization scheme is employed 
in terms of the equation of motions of the wave functions. The most popular and simplest 
one of them is the steepest descent method. The equation of motion is given by, 

(3.1) 

where is the wave function at time r and s is the Lagrange parameter. The 

damping is incorporated in Eq. 3. 1 by the following way, 

( 3 - 2 ) 

where /is the damping parameter. The main draw back of steepest descent method is 
that one can get stuck to a local minima. In another method, the damped Joannopoulos 
method [3.1], the new wave function at time r+l is calculated from those at time r-1 
and r. 

+ ro (C + )- (G + k"’ ) 

where the coefficients Vc are given by the following, 

s(h(a)-\)-{G + k\^^\G*k)e-’‘ 
e-{G + k\H^)f}*k) 
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(3.5) 


Tg =e 


-ya 


=- 


e 


h{a)-e-^^ -1 
-{G + k\H^sG + k) 


where the function h{a) is given by, 


h{a)=\ 


g-(r/2)<* 

g-(r/2)<» 


cos(iyi^) 

cosh(.^|<y^|<5) 


if co^ >Q 
if co^{0 


(3.6) 


(3.7) 


with 


6)^ =s -{G + k\H^sY^^)-r^ lA 


Though the damped Joannopoulos method is 


more efficient than the first order steepest descent method, limited usage of this method 
IS due to the fact that the process requires an additional storage for the wave function at 
time r— 1 . In order to avoid this additional storage, the William-Soler algorithm [3.2] is 
used. Here the coefficients in Eq. (2.19) are given by (with /g =0). 



s-{G+k\H,s\G->-k) 

se ' 


-{G + k\Hf,sG + k) 


(3.8) 


c-{G+k\H^\G+k) 


7c = 


s-{G + k\H^s<^ + k) 


(3.9) 


Solving the ionic equation of motion: In order to solve the atomic equation of motion, the 
Verlet algorithm [2.8] is the most common choice. The Verlet algorithm is given by, 

r(t + At) =2r(t)- r{t — At) + a(t)At^ + 0(At*) (3.10) 

where a(t) = -(1/ m)W(r(t)) is the acceleration. The truncation error here is of the order 
of At ^' . But when the velocity is calculated from this equation from. 
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v(t) = + 

2At 


(3-11) 


the error associated with this expression is of the order of A/^ rather than that of . A 
different variant of Velet algorithm, known as velocity Verlet algorithm is used to handle 
this problem. Predictor-corrector [2.9] is another commonly used algorithm for this 
purpose. In this scheme, from the positions and their time derivatives, all known at time t, 
one “predicts” the position and its time derivatives at a later time t + At, by means of 
Taylor expansion. The forces are then calculated by taking the gradient at the predicted 
positions. The difference between the “predicted” and calculated acceleration is the “error 
signal”, which is used to correct the positions. 

3.1.1 Input and Output files 

There are three main input files for this program viz. starUnp, inp.mod and the 
pseudopotential files. The formulation and generation of the pseudopotenial file using the 
flii98pp code will be discussed in the next section. The file start. inp mainly deals with the 
structure of the atoms in the concerned system while the file inp.mod deals with the 
dynamics of the atom. 

In the file inp.mod, the scheme for the generation of the wave function is 
provided. The generation of the wave function can be done by actually solving the Kohn- 
Sham equation or from file fort. 70 (in case of restarting the program or in the case of 
dynamics of the atom where the program can be started from the self-consistent electron 
densities). The step length of electronic iteration (delt) and the damping factor {gamma) 
is provided in this file. It is a common practice to choose the above two to be very large 
initially and change it to a smaller one when the self-consistency is almost reached. Apart 
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from this, various convergence criteria such as energy (epsel), kinetic energy (epsekinc) 
and force (epsfor) is to be provided in it. This program supports three different schemes 
to iterate the electronic wave functions, viz. steepest descent, damped Joannoupolos 
algorithm and William-Soler algorithm. For solving ionic equation of motion, any one of 
predictor-corrector or Verlet algorithm may be used. The exchange and correlation 
energy can be calculated by LDA or GGA approximation. 

The other input file is start. inp. The cell type (ibrav) i.e. simple cubic lattice, bcc 
lattice, fee lattice and others, and the point group symmetry (pgind) is provided in this 
file. An automatic search option of point group symmetry is also provided, which is 
useful in case of cells with unknown symmetry. The cell dimensions are to be given 
according to the choice of the cell type. The unit is atomic unit (a.u., 1 a.u. = 0.529 1 8 A. ). 
The special k-points, generated either by Chadi-Cohen [3.3] or Monkhorst-Pack [3.4] can 
be used. 

The main output file is fort. 6, which contains information about the total and 
interaction energies. A complete protocol of the run (electron energy minimization or 
molecular dynamics) is also stored in this file. Among the other output files, fort.l 
contains the instantaneous atomic positions. The binary format output files fort.71 
contains complete restart information. Other important binary format output files are 
fort.72 and fort.73, which contains the electronic charge density and total effective 
potential. Another output file fort. 21 contains the positions and velocities in binary 
format, when molecular dynamics is done. After the job is done, the following 
convergence criteria should be checked : 
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The difference between the total energy and harris energy should not be more than 
the parameter epsel provided in the file inp.mod. 

• The program should stop before the CPU time limit exceeded. 

• The program should stop before the maximum number of electronic iteration 
(nomore) or atomic movement (nstepe). 

3.2 Effect of the various parameters used in the fhi96md code 


The effect of various parameters, which control execution of calculation of total 
energy, is discussed below with a purpose of learning about the code. Here all the 
calculations are done using a single unit cell of GaAs. The plane wave cutoff energy used 
is 8 Ry and the experimental lattice constant of GaAs, 5.64191 A (10.6616 a.u.) [3.5] is 
used in all the cases. For reference, input files are attached in appendix A. 

File : Start.inp 

1. nel_exe : neljexe is the number of excess electrons present in the system. For an excess 
electron we put it as +1 while for an excess hole we put it as -1. For a neutral cell the 
total energy is -33.998000 a.u.. When an excess electron is present in the system, the total 
energy becomes -33.863005 a.u., while with an excess hole the energy is -34.059153 a.u.. 

2. ibrav : This parameter ibra\ determines the type of the lattice to be used in the 
program. As the GaAs cell is similar to zinc sulfide (interpenetrating fee) in nature, using 
an fee cell (ibrav = 2) is the best choice. For this case, we have to specify the position of 
the Ga (0,0,0) and the position of the As ( 7 , 7 , 7 ) • But, ultimately we are interested in 

incorporating vacancy defects in the cell. In that case, when we try to remove one Ga or 
one As atom from the cell, all the atoms of the given type are removed. To handle this 
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problem, the cell is chosen as simple cubic one {ibrav = 1). In this cell, 4 Ga atoms are 
placed at (0,0,0) (y,y,0) (^,0,^) (0,y,-j) while the 4 As atoms at ( 1 , 7 , 7 ) 

(l.f m) (fj iji) (i’l’i)- this case, removing one atom from a specific site, still 
keeps the remaining atoms. The transformation from fee to simple cubic structure is to be 
done in such a way that the total energy per pair should remain the same. For a cell with 
simple cubic symmetry, the total energy is -33.998000 a.u.. This energy is for 4 pairs of 
GaAs. Thus, the energy per pair is -8.4995 a.u. When the same calculation is repeated for 
a cell with fee symmetry, the energy per pair is -8.499501 a.u. The difference in energy 
of ~ 10 '^ a.u. is minor. 

3. k-point : The choice of the special k-point is a structure specific issue. The special k- 
points are generated either according to the Chadi-Cohen [3.3] mesh or according to 
Monkhorst-Pack [3.4] mesh. Since GaAs has fee symmetry, here we have used the 
special k-point for the fee structure (0.5, 0.5, 0.5). The total energy obtained is 
-33.998000 a.u.. When the same calculation is repeated for a special k-point (0.0, 0.0, 
0.0), the energy obtained is -34.002714 a.u. 

4. fold parameter : The fold parameter denotes the folding of the k points. A (3 3 3) 
folding means there is a (3>«3x3) folding of the k-point used. Accuracy of the calculation 
is better with large folding. In order to reduce the computational load, we have used a 
( 1 x 1 x 1 ) folding or no folding of the k-point in our calculation. When a (3x3x3) folding 
is used, the total energy obtained is -34.355107 a.u. (compare with -33.998000 a.u., used 
in our calculation with(l x 1 x i) folding). 

5. Cutoff energy : The effect of cutoff energy will be discussed in detail in chapter 4. 
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6. ekt : The parameter ekt denotes the temperature of the artificial Fermi smearing of the 
electrons. This is measured in eV. When this parameter is changed to 0.01 from that of its 
original value at 0.004, total energy remains the same (-33 998000 a.u ). But when this is 
changed to 1.0, total energy becomes -33.773929 a.u.. 

7. T-ion : Changing the parameter T-ton, which denotes the ionic temperature, has no 
effect on the total energy. 

8. TJnit ; The parameter TJnit, denoting the temperature of initial velocities in K, has no 
effect on the calculated total energy. 

9. Q : Q denotes the mass of the thermostat in a.u. and has no effect on the calculated 
total energy. 

File : inp.mod 

1. delt : The parameter delt determines the time step of the electronic iteration. The effect 
of the parameter delt is shown in table 3.1. 


Table 3.1 Variation of the convergence time with the parameter delt 


delt 

time taken for 
convergence (sec) 

Total Energy (a.u.) 

10.0 

38.77 

-33.997989 

20.0 

34.70 

-33.998000 

30.0 

32.32 

-33.998012 

40.0 

32.04 

-33.998012 

50.0 

28.72 

-33.998013 


2. gamma : The parameter gamma is the damping parameter used in the damped 
Joannoupolos algorithm. The effect of this parameter is shown in table 3.2. 




Table 3.2 Variation of the convergence time with the parameter g'n/n/na 


gamma 

time taken for 
convergence (sec) 

Total Energy (a.u.) 

0.5 

28.96 

-33.998013 

0.7 

29.68 

-33.998010 

0.9 

26.53 

-33.998000 

1.1 

46.26 

-33.998013 

1 3 

47.83 

-33.998014 


3. i_edyn : i_edyn denotes the scheme for the iteration of the wavefunctions. The total 
energies obtained using the steepest descent method, William Soler algorithm and 
damped Joannoupolos are -33.997984 a.u., -33.997999 a.u. and -33.998000 a.u., 
respectively. The difference in the three cases are ~10'^ a.u.. 

3.3 The fhi98pp program structure 

The program fhi98pp is needed to generate the pseudopotentials. The Fortran 
program fln98pp has two mam utilities: (i) a UNIX C-shell script psgen with related 
program fhipp and (ii) another C-shell script pswatch with associated program pslp. The 
program fliipp calculates the ionic pseudopotential for an element while pslp tests the 
transferability of the pseudopotentials thus generated. The input file (input files for the 
generation of the Ga, As, Si and C pseudopotential are attached in appendix B) is 
provided with the atomic number, number of core and valance states, atomic 
configuration and the type of exchange-correlation approximation used. Artificial atoms 
with fractional atomic number can be handled with this program. Both, Bachelet- 
Hamann-Schuter and Troullier-Martins type of pseudopotentials can be generated by this 
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program. The default cutoff radius for the pseudo wavefunction is provided with the shell 
script utility psgen. For Bachelet-Hamann-Schiiter case, it is of the outer most 

maximum of the all electron wavefunction and for Troullier-Martins, it is if some 

core states are present with the same angular momentum or otherwise. The 

default cut-off radius can be changed according to the need. Another shell script utility 
pswatch is used after the pseudopotential is generated by psgen. It performs a self- 
consistent calculation of the pseudo atom, computes the logarithmic derivatives of the all 
electron potential and pseudopotential and carries out analysis for ghost states for 
Kleinman-Bylander form of non-local pseudopotential. There are several parameters that 
control the evaluation of the logarithmic derivatives and the fully separable 
pseudopotential. 

All electron and pseudo wavefunction plots of gallium (fig.3.2-3.5), arsenic 
(fig.3.6-3.9), silicon (fig.3.10-3.13) and carbon (fig.3.14-3.17) along with their 
pseudopotentials are shown in the following section. From the plots of the wavefunction 
it is clear that all the pseudo wavefunctions thus generated are nodeless and beyond a 
certain cutoff radius, all the pseudo wavefunctions exactly matches with the 
corresponding all electron wavefunctions. As a final check, we compare the pseudo 
wavefunction eigenvalues with the corresponding all electron valance wavefunctions. All 
the results along with the cutoff radius are given in table 3.3. 
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Table.3.3 Comparison of the all electron and pseudo wavefunctions along with the cutoff 
radius for Gallium, Axsenic, Silicon and carbon 


Element 

1 value 

Cutoff radius 

Eigenvalues 



(a.u.) 

All electron 

Pseudo 

Gallium 

0 

1.2234 

-9.1703234 

(4s) 

-9.1703246 

(Is) 


1 

1.4870 

-2.7331192 

(4p) 

-2.7331193 

i2pl 

Arsenic 

0 1 

1.149 

-14.6997676 

(4s) 

-14.6997684 

(Is) 


1 

1.1776 

-5.3417404 

(4P) 

-5.3417406 

JM. 

Silicon 

0 

1.046 

-10.8791848 


-10.8791853 

(Is) 


1 

1.2714 

-4.1628555 

■SSI 

-4.1628555 

JM. 

Carbon 

0 

0.7205 

-13.6394102 

(2s) 

-13.6394110 

(Is) 


1 

0.4644 

-5.4159021 



-5.4159025 

(2p) 
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Fig.3.9. Plot for Gallium 2p pseudopotential 








Fig.3.1 1. Plot for Silicon Is pseudopotential 





Fig.3.12. Plot for Silicon pseudo 2p and all electron 3p wavefunction 
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3.15. Plot for Carbon Is pseudopotential 




Fig.3.16. Plot for Carbon pseudo 2p and all electron 2p wavefunction 
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Chapter 4 

Results and Discussion 

4.1 Optimization of the Gallium Arsenide structure 

Gallium arsenide (GaAs) has the well known zinc blande (ZnS) structure. In a 
gallium-arsenide unit cell, four gallium atoms are placed at 
(0,0,0) (^,1,0) (^,0,]-) (0,y,y) while the four arsenic atoms are placed at 

(t.T> 4 ) (Imm) (■!.■?. i) The supercell for computation is made up of 8 unit 

cells containing 32 atoms of Ga and As each. In order to optimize the stmcttire, we first 
determined the cutoff energy at which the subsequent calculations are to be done. The 
choice of the cutoff energy is a cmcial factor. Ideally, the number of plane waves used in 
this calculation should be infinite. According to energy cut-off, the number of plane 

waves are, ngwx=-^Q£l'J , where n is the volume of the supercell and is the 
6k 

energy cutoff (in Ry). The larger the number of the plane waves, the better is the 
accuracy of the results. But, a large cutoff energy requires more computational time. 
Therefore, the cutoff energy is chosen in such a way that at this energy, the variation of 
the total energy with cutoff energy does not affect the accuracy of the calculation. In 
order to establish the accuracy of our calculation, we will compare our results with a 
measured value of 0.74 eV for heat of formation of GaAs molecule. Thus, allowable 
value of energy in all minimization problems is for the condition that energy does not 
change by more than 0.01 eV per atom (molecule) (i.e. 0.001 a.u , 1 a.u. = 27.21 leV) 
upon changing the parameters in the minimization problem. 


49 



The total energy is calculated at different cutoff radii and the results are plotted to 
get an optimum value. The results are shown in table 4. 1 and the corresponding plot is 
shown in fig 4.1 

For our calculation, we have chosen the cutoff energy as 20 Ry. This choice is 
made because of following considerations. The time taken for the calculation of the total 
energy at 18, 20 and 22 Ry cutoff energy were approximately 5,7 and 12 hours, 
respectively. Though the total energy at 22 Ry is lower than that at 20 Ry cutoff, but only 
by 0.04 a.u., and calculations at 22 Ry consumed extra computational time. Also note 
that, this chosen cutoff is already higher than the equivalent calculations reported in 
literature [1-7]. Bachlet, Hamann and Schltiter [3.2] had done the calculatiori of heat of 
formation at 18 Ry cutoff and estimated it at 0.73 eV, while Chadi and others [1.7] had 
done the same calculation with a soft pseudopotential at 6 Ry cutoff and got an error of 
-0.35 eV in energy. Thus, it can be expected that for our calculation, a 20 Ry cutoff 
energy is sufficient. The input files start inp and inp.mod used for the optimization of the 
GaAs structure is given in appendix C The supercell was constmcted by assuming a 
simple cubic lattice and atomic positions of all 32 atoms of As and Ga are specified in 
start.inp. 

Once the cutoff energy is determined, the optimization for the lattice parameter is 
straightforward. Again we plot total energy against the lattice parameter. The minimum 
of the curve gives the required optimum point. The optimum value of the lattice constant 
(the lattice constant of the supercell, 2a) obtained at this cutoff is 20.98 a.u. and 
calculated total energy at this point is -277.113509 a.u. Since this energy is for 32 pairs of 
Ga and As atoms, the total energy per pair is -8.659797 a.u. (-235.64174 eV). Also, since 
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the supercell is made up of 8 actual Zinc blend GaAs unit cells, the lattice parameter for 
GaAs IS 10.49 a.u. 


Table 4. 1 Variation of total energy with the cutoff energy 


Cutoff Energy (a.u.) 

Total Energy (a.u.) 

8 

-274.712182 

10 

-275.601537 

12 

-276.320312 

14 

-276.747908 

16 

-276.953976 

18 

-277.048783 

20 

-277.110737 

22 

-277.151251 


Table 4.2 Variation of total energy with the lattice parameter 


Lattice Parameter (a.u ) 

Total Energy (a.u.) 

20 

-276.837322 

20.2 

-276.943782 

20.4 

-277.024203 

20.6 

-277.077615 

20.8 

-277.107448 

20.9 

-277.111396 

20.95 

-277.110575 

20.98 

-277.113509 

20.985 

-277.113289 

20.99 

-277.113163 

21 

-277.112678 

21.1 ■ 

-277.113194 

21.2 

-277.103183 

21.4 

-277.075901 

21.6 

-277.031815 

21.8 

-276.976025 

22 

-276.910593 









4.2 Optimization of the Arsenic structure 


The structure of arsenic can be described in three different ways [4.1]: 

(1) pnmitive rhombohedral cell with a slightly less than 60°, two atom? at ± (x,x,x) , 
where x is less than Va, 

(2) hexagonal unit cell with six atoms at [(0,0,0) ( | , | , i) (^ , f , f )]± (0,0, z) where z has 
the same value as that of x, and 

(3) face centered rhombohedral unit cell with a less than 90° and eight atoms 
at [(0,0,0)(0,Y,^)(-j,0,-j)(-j,Y,0)]±(|,-|-,'|-), where x is the same as above. 

In our calculation, the stmeture has been described according to the second definition and 
this structure of As is shown in figure 4.3. The lattice constant and the values of the 
postional parameter are taken as a = 3.7599 A (7.10515 a.u.), c = 10.5478 A (19.9323 
a.u.) and z = 0.2271 [4.1]. In actual calculation, a supercell is constructed by considering 

four unit cells of arsenic. These four cells were taken by repeating unit cells in a and b 

directions of the lattice, which mean lattice parameter taken in a and b direction is 
2x7. 10515 a.u. = 14.2103 a.u. and also the 24 atomic positions for the supercell are given 
in appendix D. Since the lattice parameters described above are not optimized ones (i.e. 
the total energy is not minimum at the corresponding values), they are optimized through 
the following procedure. The total energies were calculated with five different values of 
lattice constants a (12.7887, 13.4998, 14.2103, 14.9208 and 15.6313 a.u.), c (17.9391, 
18.9357, 19-9323, 20.9289 and 21.9255 a.u.) and z (0.2044, 0.2158, 0.2271, 0.2385 and 
0.2498) for all their possible combinations. The maximum and the minimum values of 
these parameters represent ±10% variation about the central value With the energies 
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Fig.4.3(a-d). Plot for As structure (Fig. 4.3(d) is taken from Ref [4.3]) 







calculated for all 125 sets formed by combinations of 2a, c and z, the minimum energy 
condition was obtained by a neural network program. This program trains the network for 
the given 125 sets of 2a, c and z values and corresponding energies included in appendix 
E From the allowable range of variation for the different parameters, the program then 
calculates the energy of a large number of points and predicts the point of minimum 
energy from them. The point of minimum as descnbed by the neural network program in 
this case is 2a = 15.0994 a.u., c = 20.0898 a.u. and z = 0.2748. Here the predicted z value 
IS outside the range of our input data set (z is allowed to vary from 0.2044 to 0.2498). 
This shows that total energy is a strong function of the parameter z. In order to estimate 
the accuracy of the prediction by the neural network optimization, we also graphically 
represent the energy as a function of 2a, c and z. 

Examine the plot of the total energy against the lattice parameter 2a in fig. 4.4. 
For a fix value of c, the total energy becomes independent of the value of a for 2a > 15 
a.u.. Also for 2a > 15 a.u., the variation in energy with respect to z is minimal when z is 
large. Like wise, it is clear from fig. 4.5, for 2a between 14 and 15 a.u., energies are 
insensitive to the value of c. Thus prediction of 2a = 15.0994 a.u. seems acceptable. 

From fig. 4.6, it can be seen that in the given range of value of c (17.9391 to 
21.9255 a.u.), total energy is insensitive to the variation in the value of z when z is large. 
Thus the z dependency of the total energy may be eliminated at a higher value of z. The 
next curve, fig. 4.7, indicates that, for 2a between 14.2 and 14.9 a.u., total energy does 
not change significantly with an increase in the value of c. Since these calculations are 
done at z = 0.2498, it is expected that c dependency of the total energy will be less at a 
still higher z value. 
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From plots in fig 4.6 and fig. 4.7, c — 20.0898 a.u. does not appear unreasonable 
as a first estimate 

However, a more serious concern is with the prediction of neural network minima 
for z. But, it IS clear from fig. 4.8 and fig. 4.9, only for large values of z, a minimum 
energy condition is reached. Nonetheless, because the predicted value of z lies outside the 
range of input data, further investigation is necessary 

For this purpose the energy calculation is done in the neighborhood of the 
approximate minimum point, by varying one parameter at a time. If the calculated 
energies at the neighboring points are higher than the predicted point, then the predicted 
value from neural network is accepted. If the energy decreases in any given direction, 
then further calculations are performed in that direction, till either the energy increases or 
the variation of the energy with respect to the change in the parameter is 10'^ a.u. or less. 


Table 4.3 Calculated total energies at the neighboring points of the approximate 


minimum as predicted by the neural network 


Serial 

No. 

2a (a.u.) 

c (a.u.) 

z 

Total Energy (a.u.) 

1. 

15.0994 

20.0898 

0.2748 

-152.334489 

2. 

14.9994 

20.0898 

0.2748 

-152.324711 

3. 

15.1994 

20.0898 

0.2748 

-152.341461 

4 

15.0994 

19.9898 


-152.334584 

5. 

15.0994 

20.1898 


-152.334490 


15.0994 

20.0898 


-152 394685 

7. 

15.0994 

20.0898 

0.2848 

-152.166115 

8 

15.0994 

20.0898 

0.2548 

-152.392038 
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Fig.4.5. Plot for variation of the total energy with 2a for a fixed z and a family of 
values of c 
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Fig.4.6. Plot for variation of the total energy with c for fixed value of 2a and a 
family of values of z 
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Fig.4.7. 


Plot for variation of the total energy with c for a fixed value of z and a 
family of values of 2a 





Fig.4.8. Plot for variation of the total energy with z for a fixed value of c and a 
family of values of 2a 








From the data in table 4.3, it is clear that the total energy value decreases with the 
increase in the value of parameter 2a. The results show that if the parameter a is 
increased, the total energy decreases. But this decrease in energy is ~10'^ a.u./atom. A 
variation of IQ-' a.u. in total energy, which is equivalent to lO'^ eV, is insignificant in our 
calculation. No such optimization of the parameter c is required, as when this parameter 
IS changed, total energy of the arsenic remains either constant or increases from that at 
the neural network minimum (see rows 1, 4 and 5 in table 4.3). Also it is clear from the 
table 4.3 that the total energy is a strong function of z and needs special attention. We 
find a more suitable minimum for z at 0.2648 (see rows 1, 6, 7 and 8 in table 4.3). Thus, 
from this consideration, the point of minimum is fixed at 2a = 15.0994 a.u., c'= 20.0898 
a.u. and z = 0.2648. We once again repeat the procedure adopted to generate table 4.3 to 
examine if the sensitivity to a and c has changed at this new value of z. The results are 
shown in table 4.4. Clearly, 2a = 15.0994 a.u., c = 20 0898 a.u. and z = 0.2648 can be 
accepted as the lattice parameters. Thus the minimum energy for this 24 atom arsenic 
supercell is -152.394685 a.u. (i.e. -172.783836 eV/atom). 

Table 4.4 Total energies around the point of minimum energy 


Serial 

No. 

2a (a.u.) 

c (a u.) 

z 

Total Energy (a.u.) 

1. 

15.0994 

20.0898 

0.2648 

-152.394685 

2. 

15.0994 

20.0898 

0.2848 

-152.166115 

3. 

15.0994 

20.0898 

0.2548 

-152.392038 

4. 

15.1994 

20.0898 

0.2648 

-152.392158 

5. 

14.9994 

20.0898 

0.2648 

-152.396238 

6. 

15.0994 

20.1898 

0.2648 

-152.391698 

7. 

15.0994 

19.9898 

0.2648 

-152.398110 
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4.3 Optimization of the Gallium structure 


Optimization of the gallium is done in the same way as that of arsenic. Gallium, at 
room temperature, has orthorhombic symmetry with two of the axes of nearly equal 

length. The eight gallium atoms are placed at [(0,0,0) (-,-,0)]±[(0,y,z) (-, y, --z)]. 

This structure is plotted in tig. 4.10. At room temperature the lattice parameters are a = 
4.5192 A (8.54 a.u.), b = 7.6586A (14.4726 a.u.) and c = 4.5258 A (8.5525 a.u.) while the 
offset parameters y and z are 0.1539 and 0.0798 [4.1], respectively. Here the supercell is 
constructed by considering four gallium unit cells, and therefore, the supercell contains 
32 gallium atoms. The lattice parameters of the supercell for computation taken in a and 

b direction are 2x8.54 = 17.08 a.u and 2x14.4726 = 28.9452 a.u. respectively. In order 
to optimize the structure, the above said structural parameters (a, b, c, y and z) are varied 
by ±10% about the room temperature value given above. Total energies are calculated for 
all the combinations of three different values of 2a (15.372, 17.08 and 18.788 a.u.), 2b 
(26.0507, 28.9452 and 31.8397 a.u.) and c (7 6973, 8.5525 and 9.4078 a.u.) and two 
different values of _!■ (0.1385 and 0.1693) and z (0.0718 and 0.0878). The input file for 
fhi96md, which also provides the 32 atom positions, is attached in appendix F. The 
calculated value of total energies, for vanous combinations of 2a, 2b, c, y and z, as 
mentioned above, are given in appendix G. These 108 set of data containing 2a, 2b, c, y, 
z and total energy are then fed to a neural network program to get the point of minimum 
energy for a given set of lattice parameters. The point of minimum energy is predicted by 
the neural network at 2a = 16.514 a.u., 2b = 28.8253 a.u., c = 8.0075 a.u., y = 0.1586 and 
z = 0.0566. 
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Total energies are calculated at the neighboring points of this approximate 
minimum point by varying one parameter at a time. The results are shown in table 4.5. 


Table 4.5 Calculated total energies at the neighboring points of the approximate 
minimum predicted by the neural network 


a (a.u.) 

b (a.u.) 

c (a.u.) 

y 

z 

Total Energy 
(a.u.) 

16.5140 

28.8253 

8.0075 

0.1586 

0.0566 

-73 398999 

16.6140 

28.8253 


0.1586 

0.0566 

-73.390369 

16.4140 

28.8253 

limM 

0.1586 

0.0566 

-73.407124 

16.5140 

28.9253 


0.1586 

0.0566 

-73.396928 

16.5140 

28.7253 


0.1586 

0.0566 

-73.400781 

16.5140 

28.8253 


0.1586 

0.0566 

-73.379983 

16.5140 

28.8253 

7.9075 

0.1586 

0.0566 

-73.415477 

16.5140 

28.8253 

8.0075 

0.1686 

0.0566 

-73.291836 

16.5140 

28.8253 


0.1486 


-73.360814 

16.5140 

28.8253 

8.0075 



-73 382359 

16.5140 

28.8253 

8.0075 

0.1586 

0.0466 

-73.411612 


Clearly, with both increase and decrease in the value of y about the predicted minimum, 
the cell energy decreases Thus the predicted minimum for y is acceptable. Whereas, in 
the case of parameters b, the energy increases with an increase in the value of b. 
However, for a reasonable change in the value of b, the change in energy is of the order 
of 10'^ a.u. Hence, we accept the predicted minimum for b. However, a change in cell 
energy as a function of a, c and z is significant. Thus, a better estimate for a, c and z is 
searched further. At this stage, we vary only one of the parameters a, c and z at a time. 
To get an optimum value of c, total energies are calculated at different values of c 
keeping the other parameters constant at predicted minimum of the neural network The 
corresponding values are shown in table 4.6. The results are also shown in graphical form 
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m fig 4.11. It shows a distinct minimum at c = 7.6075 a.u., which is now taken as the 
new optimum value. In order to get the optimum value of z, the same procedure is 
applied. Total energies are calculated at different points by decreasing the value of z and 
keeping the other parameters constant at value shown in fig. 4.12. The results are given in 
table 4.7 and the corresponding graph is shown in fig. 4.12. In this case, instead of a well 
defined minimum, we get a curve that decreases with the decrease in the value of z. Since 
a variation of 10 ' a.u. in total energy is allowed in our calculation, we chose the point z = 
0.0366 as the optimum point. Similarly, the total energy value also decreases with a 
decrease in parameter a. As no distinct minimum was obtained in this case also, the 
optimum point is chosen at 2a = 16.1 14 a.u. The variation of total energy with the lattice 
parameter a is given in table 4.8 and is plotted in fig. 4.13. Thus, in subsequent 
calculations, the following optimum values are taken: 2a = 16.114 a.u., 2b = 28.8253 a.u., 
c = 7.6075 a.u., y = 0. 1 586 and z = 0.0366. 


'fable 4.6 Variation of total energy with parameter c 


c (a.u.) 

Total Energy (a.u ) 

8.1075 

-73.379983 

8.0075 

-73.398999 

7.9075 

-73.415477 

7.8075 

-73.427118 

7.7075 

-73.433446 

7.6075 

-73.433953 

7.5075 

-73.428174 
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Fig.4.1 1. Plot for total energy of gallium as a function of parameter c 



Table 4.7 Variation of total energy with parameter z 


z 

Total Energy (a.u.) 

0.0798 

-73.381665 

0.0566 

-73.433953 

0.0466 

-73.448832 

0.0366 

-73.459965 

0.0266 

-73.467873 

0.0166 

-73.473001 

0.0066 

-73.475657 


Table 4.8 Variation of total energy with parameter 


a (a.u.) 

Total Energy (a.u.) 

16.514 

-73.459965 

16.414 

-73.469243 

16.314 

-73.477905 

16.214 

-73.485893 

16.114 

-73.493255 

16.014 

-73.499882 

15.914 

-73.505717 

15.814 

-73.510774 
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Fig.4.12. Plot for total energy of gallium as a function of parameter z 
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Fig.4. 1 3 . Plot for total energy of gallium as a function of parameter a 




Finally, we examine the sensitivity of energy around this new point of minimum. These 
results are shown in table 4.9. 


Table 4.9 Total energies around the point of minimum energy 


2a (a.u.) 

2b (a.u.) 

c (a.u.) 

y 

Z 

Total energy 
(a.u.) 

16.1140 

28.5253 

7.6075 

0.1586 

0.0366 

-73.495423 

16.0140 

28.5253 

7.6075 

0.1586 

0.0366 

-73.501631 

16.2140 

28.5253 

7.6075 

0.1586 

0.0366 

-73.488497 

16.1140 

28.4253 

7.6075 

0.1586 

0.0366 

-73.487144 

16.1140 j 

28.6253 

7.6075 

0.1586 

0.0366 

-73.495052 

16.1140 

28.5253 

7.5075 

0.1586 

0.0366 

-73.491861 

16.1140 ! 

28.5253 

7.7075 

0.1586 

0.0366 

-73.492762 

16.1140 

28.5253 

7.6075 

0.1486 

0.0366 

-73.461380 

16.1140 

28.5253 

7.6075 

0.1686 

0.0366 

-73.349660 

16.1140 

28.5253 

7.6075 

0.1586 

0.0266 

-73.503096 

16.1140 

28.5253 

7.6075 

0.1586 

0.0466 

-73.484418 


It is evident from the results in table 4.9 that the new values of cell parameters are either 
at a minima, or energy does not change significantly when any of the parameters are 
changed. The calculated total energy at the point of minimum is -73.495423 a.u.. Since 
this energy is calculated for a supercell containing 32 Ga atoms, the total energy per atom 
is -2.296731 a.u. (-62.496374 eV). And the final optimized lattice parameters obtained 
are a = 8.57 a.u., b = 14.26265 a.u., c = 7.6075 a.u., y = 0.1586 and z = 0.0366. 
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4.4 Heat of formation of GaAs 


The bulk energies of Ga, As and GaAs are calculated, as mentioned in Section 4. 1-4.3, 
using optimized lattice parameters. To calculate the heat of formation of GaAs, the bulk 
energies of Ga and As are subtracted from the total bulk energy of GaAs (per pair). The 
calculated heat of formation of GaAs is compared with a set of published results [1.7], 
where various pseudopotentials were used in the calculation. The table 4.10 contains 
those total energy and heat of formation values along with the experimental heat of 
formation of GaAs [4.2] and our calculated result. The results show that the hard 
pseudopotentials generated according to the Bachelet, Hamann and Schlutef [2.11] and 
Kerker [1.8] scheme (0.72 and 0.73 eV respectively) matches well with the experimental 
result, while the soft pseudo potential gives the heat of formation as 1.08 eV. When these 
calculations are repeated for both hard and soft pseudopotentials generated according to 
the Kerker scheme at 6 Ry cutoff energy (modified Kerker scheme as applied by Chadi 
and others [1.7]), the heat of formation obtained to be 1.76 eV and 1.18 eV respectively. 
Clearly, the heat of formation for the hard pseudopotential is ~1 eV off from that of at 12 
Ry. On the other hand, the soft pseudopotential result is off by 0.1 eV. In order to 
compensate this error at 6 Ry, this soft pseudopotential results are corrected according to 
the following considerations. The soft pseudopotential heat of formation (1.08 eV) is 
-0.35 eV off from the expenmental result. To correct this, the Ga and As heat reservoirs 
are lowered by an amount of 0.175 eV (=0.35/2 eV), while keeping the energy of GaAs 
fixed. The other error of 0.1 eV between the 12Ry and 6 Ry cutoff energy is equally 
distributed between Ga and As. This corrected result is also included in the table 4.10. 
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The calculated bulk energies of Ga, As and GaAs, in our case, are -62.496374, - 
172.783836 and -235.64174, respectively. The heat of formation obtained is 0 36 eV, 


Table 4.10 Bulk energies of Ga, As and GaAs and the heat of formation of GaAs. 



Cutoff 

energy 

(Ry) 

Bulk Energy (eV) 

Heat of 
formation 
(eV) 

GaAs 

Ga 

As 

Measured* 





0.74 

Hard Ga, As'^ 

18 

-235.379 

-61.364 

-173.294 

0.72 

Hard Ga, As^ 

12 

-235.832 

-62.613 

-172.490 

0.73 

6 

-230.825 

-60.988 

-168.075 

1.76 

SoflGa,As^ 

12 

-235.512 

-62.364 

-174.070 

1.08 

6 

-236.168 

■ -62.161 

-172.823 

1.18 

Soft Ga, As 

12 

-237.512 

-62.538 

-174.245 


Soft Ga, As^ 

6 

-236.168 

-62.388 

-173.051 

0.73 

Present 

calculation 

20 

-235.64174 

-62 496374 

-172.783836 

0.36 


Experimental result [4 2] 

^ Using Bachelet, Hamann and Schluter pseudopotential at 18 Ry cutoff energy [1 7] 
^ Potential generated according to the Kerker scheme [1 7] 

* Conected heat of formation at 6 Ry 


which is off by -0.4 eV from the experimental value. On the other hand, the calculated 
bulk energies of these materials are close to the corresponding published results. We have 
also calculated the heat formation for the optimized lattice parameter of GaAs (10.49 
a.u.) and experimental lattice parameters of Ga (a = 8.54 a.u., b = 14.4726 a.u., c — 
8.5525 a.u., y = 0.1539 and z = 0.0798) and As (a = 7.1051 a.u., c = 19.9323 a.u. and z = 
0.2271). Then the bulk energies obtained are -235.64174 eV/molecule, -62.204163 
eV/atom and -172.58352 eV/atom. In this case the heat of formation obtained is 0.85 eV, 
which is very close to the experimental value. Thus the error in our calculation may be 
attributed to our adopted lattice parameter optimization procedure of Ga and As. Chadi 
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and others [1.7] had optimized the Ga structure in a different way. They had kept the ratio 
a.b.c fixed and then minimized the total energy with respect to the volume and internal 
parameters y and z. In this way they had reduced the number of vanable from 5 (a, b, c, y 
and z) to 3 (a/b/c, y and z). Our optimization procedure is more intensive one where all 
the 5 variables are allowed to change. 


4.5 Defect formation energy in GaAs 


We have investigated a Ga vacancy defect in this work. A defect is introduced in 
the perfect GaAs crystal by physically removing one atom from its lattice. Here the 
supercell for computation is constructed from 8 unit cells of GaAs. All the subsequent 
calculations are carried out at the optimized lattice parameter of 10.49 a.u. and at a plane 
wave cutoff energy of 20 Ry. The special k-point chosen for the calculation is (0.5, 0.5, 
0.5). 

In our calculation, the supercell is periodic in nature and thus, every defect 
associated with a single supercell gets repeated. When a defect is created in a single unit 
cell of GaAs by removing one Ga atom, the defect concentration becomes too high (i.e. 
we are removing one Ga atom from a GaAs cell containing 4 Ga and 4 As atoms). But we 
are interested in one vacancy in an infinite lattice. Therefore, all our defect related 
calculations are based on big supercell containing 32 Ga and 32 As atoms. Removing a 
single atom from the interior of this big supercell causes repetition of vacancy at a larger 
distance. 

In principle, the total energy of a cell containing a defect should be independent 
of the position of the defect. In order to check this, first one Ga vacancy defect is created 
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at (0.5, 0.5, 0.5) and the total energy is calculated for this supercell. The total energy 
calculated is -274.725388 a.u.. The same calculation is repeated for a Ga vacancy defect 
at (0.75, 0.75, 0.5) and in that case the energy obtained is -274.725391 a.u.. Since the 
difference between these two results are ~10 ^ a.u., this error can be neglected. Hence, in 
all our calculations, Ga vacancy defects are created at (0.5 0.5 0.5). 

A vacancy defect, at a Ga site in neutral charge state is denoted by . When an 
extra electron is present at the vacancy, created by the removal of the Ga atom, it is 
denoted by and so on. Calculated total energies of some Ga vacancy defect cells in 

different charge states along with the corresponding perfect cell are shown in the table 
4.11. 

Table 4. 1 1 Calculated total energies of a defect supercell at different charge states along 

with the perfect cell 


Cell type 

Total energy (a.u.) 

Perfect cell 

-277.113509 

V" 

-274.725388 

V' 

'^Ga 

-274.621469 

''Ga 

-274.511415 


The defect formation energy is then calculated according to the Eq. 2 41. 

= + (2-41) 

where, 

^'d - + ^As)^Gas(bulk) ~ As {bulk)) (2.42) 
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(f^Ca(bulk) f^As(bua)) (2-43) 

Ed is the total energy of the defect cell and for a cell containing a Ga defect, is 31 

22. Also as described in chapter 1, the chemical potential of the bulk crystal 

may be equated to the bulk energy per atom (molecule) of the corresponding system. The 
bulk chemical potentials, as calculated in section 4.4, of GaAs, Ga and As 
are- 235.64174 eV, -62.496374 eV and -172.783836 eV. Accordingly, our calculated 
results of the formation energies are shown in table 4.12. Our calculated results are then 
compared with the results obtained by Northrup and Zhang [2.16]. 

Table 4.12 Calculated defect formation energies of GaAs containing a Ga vacancy defect 
at different charge states along with the ionization energies 


Defect 

type 

Q(eV) 

Ionization energy (eV) 


From Reference 
2.16 

Present calculation 

From 

Reference 

2.16 

Present 

calculation 

yO 

4.55 + Ap/2 

2.306 + Ap/ 2 



V’ 

4.74-p.e+Ap/2 

5.1337-(pe + Ev) + Ap/2 

0.19 

2.8277 -Ev 

'^Ga 

4.94-2pe+Ap/2 

8.128 -2(Pe+Ev) + Ap/2 

0.2 

2.9943 -Ey 


Here, pe can be considered as fermi energy measured from the valence band edge. 
Therefore, in order to compare our results with that of Northrup and Zhang [2.16], we 
need to estimate Ey. Furthermore, since calculations of ionization energy involves 
difference in formation energies of same defect in different charge states, better match 
with those in [2.16] can be expected. 
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The ionization energy is calculated according to the following consideration. A 
defect changes its ionization state in order to minimize the formation energy. For 
example, in our calculations, we find for fermi energy (pe) 2.8277-Ev (eV), the formation 

energy of and Vga is the same. For pe less than this value, is stable and for pe 

greater, Vq^ is more stable. Thus, the acceptor level (0/1) due to Ga vacancy lies 
2.8277 - Ev (eV) above Ev. 

Unfortunately, however the Ey energy level for calculation of total energy of a 
perfect and a defect cell are not the same. And as a result, a band offset correction to this 
valance band top is needed. 

In order to apply these corrections, band stmcture calculation is done for the T 
point (k =0, 0, 0) and the eigenvalue of the highest occupied state is noted. This is the 
energy of the top of the valance band and in our case it is 2.884 eV. Again, there may be 
an overall offset of the electrostatic potential that is different in the bulk cell and in the 
defect cell. To correct for this, the difference in the eigenvalues of the lowest occupied 
states of these two states are used as a correction. Thus corrected Ey is given by, 

Evuorrece^) = Ey - (s ^ ^ where s and urc eigenvalues of the lowest 
occupied states of the perfect and the defect cell calculated at k=(0.5, 0.5, 0.5). The 
calculated value of £ perfect ^defect (^Ga^ -9.505 and -9.694 eV, respectively. 

Thus, corrected Ey is at 2.695 eV. Accordingly, the ionization levels (0/-1) and (-1/-2) for 
Ga vacancy are at 0.1227 and 0.2893 eV above the valence band edge. 

Poykko and others [4.4] had discussed about this correction in detail. According 
to them, the difference in the energy level of the perfect supercell and the defect supercell 
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IS due to the finite size of the cells. The reference energy level value (which is the the 
valence-band maximum, VBM) is then corrected firom the difference of the average 
potentials of the defect supercell and the perfect supercell, i.e. 

Ey ^ V {bulk) average{bulklike) ^ average{bulk) (^*^) 

For a very large supercell, the last two terms in Eq. (4.1) cancels out. However, using 
such a big supercell for this type of supercell calculation is impossible. 

Lakes and others [4.5] had also discussed the effect of the presence of a defect on 
the reference energy level of total energy calculation. They had used the “model solid” 
theory as proposed by Van de Walle and Martin [4.6] to calculate the average potential of 
a system and thus to predict the shift in energy due to the incorporation of the defect. 
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Chapter 5 

Conclusions 

In the present work we have studied the defect formation energies of gallium 
vacancies in GaAs using the computer code fhi96md. In the process, we have 
generated first principle, fully separable norm conserving pseudopotentials using the 
associated code fhi98pp 

The pseudopotentials for Ga, As, Si and C are generated according to Bachlet- 
Hamann-Schliiter scheme. To ensure the accuracy of the pseudopotentials, all electron 
valance eigenvaules were compared with the corresponding pseudo wavefiinctions. 

The effects of various parameters appearing in the fhi96md code were 
examined. In order to optimise the GaAs structure, the effect of the two most 
important parameters, cutoff energy and lattice parameter, were discussed in detail 
On the basis of the total energy minimization, the optimised values of cutoff radius 
and lattice parameter was chosen as 20 Ry and 10.49 a.u., respectively. These two 
parameters were kept fixed for the subsequent calculations. 

To optimise the structures of Ga and As, a neural network program was used. 
This program gives an approximate point of minimum for the total energy. To check 
the accuracy of the neural network program, some of the inputs of the neural network 
program used for the optimization of Ga, were plotted. These plots show the same 
trend. Finally, total energy calculations were done at the neighbour-hood of the 
approximate point of minimum. Again, on the basis of the total energy minimization, 
the optimised lattice parameters for these two stmetures were chosen. The optimised 
lattice parameters of Ga is chosen as a = 8.57 a.u., b = 14.26265 a.u , c = 7.6075 a u , 
y = 0.1586 and z = 0.0366. For Ga, the total energy per atom, which is also the 
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chemical potential, is given by, -62.496374 eV. Similarly, for As, the optimised lattice 
parameters are given by, a = 7.5497 a.u., c = 20.0898 a.u. and z = 0 2648. The 
chemical potential for As is estimated as -172.783836 eV 

The heat of fomiation for these optimised stmctures of GaAs, Ga and As are 
then calculated and compared with the results published in the literature. Our result 
(0.36 eV) is less than the experimental measure (0.74 eV). The discrepancy is mainly 
due to the optimization procedure we have adopted to optimize the stmcture of Ga 
and As. 

Defects were introduced to the perfect GaAs crystal and their total energies 
were calculated The calculations were repeated for different charge states of a Ga 
vacancy defect. To align the reference energy levels of the perfect and the defected 
cells, band offset correction was applied. Finally the ionization energy of the various 
defects were calculated and compared with the literature values. 
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Appendix A 

Input file used for the optimization of GaAs 

#! /bin/csh -xvf 


####### define home directory of fhi96md ############ 
set FHI96MD = "~/fln96md" 


############## set directories : 


set PSEUDO = "${FHI96MD}/pseudo" 
set SOURCES = "S{FHI96MD}/src" 
set WORK - "S |FHI96MD}/GaAs/work’’ 


########### change to the working space ############# 
cd SWORK 


#### move pscudopotentials to the working space #### 

cp $PSEUDO/ga_aa.cpi fort. 11 
cp $PSEUDO/as_aa.cpi fort. 12 


# make input file for cp-program: inp.mod 

# (controlling main features of the fhi96md) 

mumnmmwtm 


# 

# 


cat > inp.mod « cnd_B 
-1 100 1000000 
100 1 

20.0 0.9 

4.0 0.2 0.0001 
400 2 

0.0 1.0 

22 

0 -false. 

.F. 0.001 .F. 0.002 
.false, .false, .false. 1 800 
.false. 

0.0001 0.0005 0.1 
0.001 0.001 3 

1 

end B 


: nbeg iprint timequeue 
: nomore nomore_init 
; delt gamma 

: delt2 gamma2 eps_chg_dlt 
; delt_ion nOrder 
: pfft_store mesh_accuracy 
: idyn i_edyn 
; i_xc t_postc 
: trane ampre tranp amprp 
: tfor tdyn tsdp nstepe 
: tdipol 

: epsel epsfor epsekinc 
: force_eps max_no_force 
: init basis 
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J f. J t, ^ t* •J-J- 4^4 4 4 -r / iij 

rrn 7^ tt ^r fr tt *i i * ^-jx 


# create input file : inp.mi 

# (controlling the atomic structure and the input) 

# and parameter.h (used for compiling of £hi96md) 

# by the help of the program start. f 

# make inp for start program 


# 

# 

# 

# 

# 


cat > start, inp « end_ 

_A ##################### 


2 

: number of species 


0 

: excess electrons 


5 

: number of empty states 


1 0 

: ibrav pgind 


10.6616 0.0 0.0 0 00 

: celldm 


1 

: number of k-points 


0.5 0.5 0.5 1.0 

: k-point coordintes, weight 


1 1 1 

: fold parameter 


.false 

: k-point coordinates relative or absolut? 


8 4.0 

; Ecut [Ry], Ecuti [Ry] 


0.004 .true. .f. 

: ekt tmetal tdegen 


.true, .false. 1 

; tmold tband nrho 


5 2 1234 

: npos nthm nseed 


873.0 1400.0 le8 1 

: T_ion T_init Q nfi_rescale 


.t. .true. 

: tpsmesh coordwave 


4 3 'Gallium' 

: number of atoms, zv, name 


1.0 69.72 0.7 3 3 

: gauss radius, mass, damping,l_max,l_loc 


.t. .t. .f 

: t_init_basis 


0.0 0.0 0.0 .t. 

: tauO tford 


0.5 0.5 0.0 .t. 

: tauO tford 


0.5 0.0 0.5 .t. 

: tauO tford 


0.0 0.5 0.5 .t. 

: tauO tford 


4 5 'Arsenic' 

: number of atoms, zv, name 


1.0 74.92 0.7 3 3 

: gauss radius, mass, damping, l_max,l_loc 


■t. .t. .f 

: t_init_basis 


0.25 0.25 0.25 .t. 

; tauO tford 


0.75 0.25 0.75 .t. 

: tauO tford 


0.75 0.75 0.25 .t. 

: tauO tford 


0.25 0.75 0.75 .t. 

: tauO tford 


end_A 



# run start program and build up parameter.h 

# 

# 

# and inp. ini 


# 

# -in principle one could create inp.ini and 

# 

# parameter.h by hand, but start gives a 

# 

# consistent and optimized input for fhi96md 

# 

■# 



pushd ${SOURCES}/start 
make 

if (Sstatus) then 

echo 'makes.start: compilation failed' 
exit 
endif 
popd 

cp ${SOURCES}/start.x . 
start.x 1 tee start. out 


# 



compile flii96md 


# 


cp $ {SOURCES }/fhi96md/parameter.h parameter.old 
set changes = 'cmp paramcter.h parameter.old' 
pushd $ {SOURCES }/fhi96md 


if( "Schanges" != ) then 

cp SWORK/parameter.h . 
make clean 
endif 


make netlib 


if (Sstatus) then 

echo 'makcs.flii96md: compilation failed' 
exit 
endif 
popd 

mv $ {SOURCES }/fhi96md.x . 




# 


run fhi96md 




T rrn nt-^irnr 


# 


fhi96md.x 



Appendix B 


Input file used for the generation of the Ga pseudopotential 

31.00 6 2 8 0.00 : z nc iiv iexc mlc 

1 0 2.00 • n 1 f 

2 0 2.00 
2 1 600 
3 0 2.00 
3 1 6.00 
3 2 10.00 
4 0 2.00 
4 1 1.00 

2 h : Imiix s_pp_def 

0 1 .25 0.00 h : It ret et sjpp_type 
2 2.25 0.00 h 


Input file used for the generation of the As pseudopotential 

33.00 6 2 8 0.00 ; z nc nv iexc mlc 

1 0 2.00 ; n 1 f 

2 0 2.00 
2 1 6.00 
3 0 2.00 
3 1 6.00 

3 2 10.00 

4 0 2.00 
4 1 3.00 

2 h : Imax s_pp_def 

0 1.15 0.00 h : It ret ct sjppjype 
2 2.05 0.00 h 


Input file used for the generation of the Si pseudopotential 


14.00 3 2 8 0.00 : z nc nv iexc mlc 

1 0 2.00 ; n 1 f 

2 0 2.00 
2 1 6.00 
3 0 2.00 
3 1 2.00 

2 h : Imax sj)p_def 
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Input file used for the generation of the C 


6. 

00 

1 

2 8 0.00 : z 

nc nv iexc mlc 


1 

0 

2.00 :n 1 

f 


2 

0 

2.00 



2 

1 

2 00 


2 

h 


: !max 

sjip_def 

0 

0. 

00 

0 00 h : It 

ret et s_pp_type 


pseudopotential 
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Appendix C 


Files startAnp and inp.mod used in the total energy calculation of GaAs 


start inp 



2 


: number of species 

0 


: excess electrons 

5 


: number of empty states 

1 0 


: ibrav pgind 

20.98 0 0 0.0 0 0 0 

: celldm 

1 


: number of k-points 

05 05 0.5 1 0 

: k-point coordmtes, weight 

1 1 1 


: fold parameter 

■false. 


: k-point coordinates relative or absolut? 

20 4.0 


: Ecut [Ry], Ecuti [Ry] 

0.004 .true. .f. 


: ekt tmetal tdegen 

.true, .false. 1 


: tmold tband nrho 

5 2 1234 


: npos nthm nseed 

873.0 1400.0 le8 1 

: T_ion T_init Q nfi_rescale 

.t. .true. 


: tpsmesh coordwave 

32 3 ’Gallium’ 


: number of atoms, zv, name 

1.0 69.72 0.7 3 3 

: gauss radius, mass, damping, l_max,l_loc 

.t. .t. .f 


• t_init_basis 

0.0 0.0 0.0 

.t. 

: tauO tford 

0.25 0.25 0.0 

.t. 

: tauO tford 

0.25 0.0 0.25 

.t. 

: tauO tford 

0.0 0.25 0.25 

.t. 

: tauO tford 

0.5 0.0 0.0 

•t. 

: tauO tford 

0.75 0 25 0.0 

.t. 

: tauO tford 

0.75 0.0 0.25 

.t. 

: tauO tford 

0.5 0.25 0.25 

.t. 

: tauO tford 

0.0 0.0 0.5 

.t. 

: tauO tford 

0.25 0.25 0.5 

.t. 

: tauO tford 

0.0 0.25 0.75 

.t. 

: tauO tford 

0.25 0.0 0.75 

.t. 

: tauO tford 

0.0 0.5 0.0 

t. 

: tauO tford 

0.25 0.75 0.0 

.t. 

: tauO tford 

0.25 0.5 0.25 

.t. 

: tauO tford 

0.0 0.75 0.25 

.t. 

: tauO tford 

0.5 0.5 0.0 

.t. 

; tauO tford 

0.75 0.75 0.0 

.t. 

: tauO tford 

0.75 0.5 0.25 

.t. 

: tauO tford 

0.5 0.75 0.25 

.t. 

: tauO tford 

0.5 0.0 0.5 

.t. 

; tauO tford 
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0.75 0 25 0.5 .t. : tauO tford 

0.75 0.0 0 75 .t. : tauO tford 

0.5 0.25 0.75 .t. : tauO tford 

0.0 0.5 0.5 .t. : tauO tford 

0.25 0.75 0.5 .t. ; tauO tford 

0.25 0.5 0.75 .t. ; tauO tford 

0.0 0.75 0.75 .t • tauO tford 

0.5 0.5 0.5 t. : tauO tford 

0.75 0.75 0.5 .t : tauO tford 

0.75 0.5 0.75 .t. : tauO tford 

0.5 0.75 0.75 .t. : tauO tford 

32 5 'Arsenic' : number of atoms, zv, name 

1 .0 74.92 0.7 3 3 : gauss radius, mass, damping, l_max,l_loc 

.t. .t. .f : t_init_basis 

0.125 0.125 0.125 t. : tauO tford 

0.375 0.375 0.125 .t. : tauO tford 

0.375 0.125 0.375 .t : tauO tford 

0.125 0.375 0.375 .t . tauO tford 

0.625 0.125 0.125 .t. : tauO tford 

0.875 0.375 0.125 .t. : tauO tford 

0.875 0.125 0.375 .t : tauO tford 

0.625 0.375 0.375 t. : tauO tford 

0.125 0.125 0.625 .t. : tauO tford 

0.375 0.375 0.625 .t. : tauO tford 

0.125 0.375 0.875 .t ; tauO tford 

0.375 0.125 0.875 .t. : tauO tford 

0.125 0.625 0.125 .t : tauO tford 

0.375 0.875 0.125 .t. : tauO tford 

0.375 0.625 0.375 .t. : tauO tford 

0.125 0.875 0.375 .t : tauO tford 

0.625 0.625 0.125 .t. ; tauO tford 

0.875 0.875 0.125 .t. : tauO tford 

0.875 0.625 0.375 .t : tauO tford 

0.625 0.875 0.375 t. : tauO tford 

0.625 0.125 0.625 .t. : tauO tford 

0.875 0.375 0.625 .t. : tauO tford 

0.875 0.125 0.875 .t. : tauO tford 

0.625 0.375 0.875 .t. : tauO tford 

0.125 0.625 0.625 .t : tauO tford 

0.375 0.875 0.625 t. : tauO tford 

0.375 0.625 0.875 .t. ; tauO tford 

0.125 0.875 0.875 .t. : tauO tford 

0.625 0.625 0.625 .t. : tauO tford 

0.875 0.875 0.625 .t ; tauO tford 

0.875 0.625 0.875 .t. : tauO tford 

0.625 0.875 0.875 .t. : tauO tford 
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inp mod 



-1 100 lOOOOOO 
lOUO 1 
20 0 0.0 
4.0 0.2 0.0001 
400 2 
0.0 1.0 
2 2 

0 false. 

F 0.001 .F. 0.002 
.false, .false, .false. 1800 
false. 

0 0001 0.0005 0.1 

0 001 0 001 3 

1 


: nbeg iprint timequeue 
: nomore nomore_init 
; delt gamma 

: delt2 gamma2 eps_chg_dlt 
: delt_ion nOrder 
; pfft_store mesh_accuracy 
: idyn i_edyn 
: i_xc t_postc 
; trane ampre tranp amprp 
: tfor tdyn tsdp nstepe 
; tdipol 

; epsel epsfor epsekinc 
: force_eps max_no_force 
: init basis 
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Appendix D 

Input file used for the o[)tiinisation of As 


;.‘l bin esh wf 

dctnic lu'iiic diiccturv’ of fln96md ############ 
SCI FHI'HAli) " Ihioomd" 

set directories 


set PSF.l 'DO ”S lFH!%MDf pseudo" 
sctSOURCl-S "S!FHl%MI>f'src" 
set Wt )RK '■$ J FH P)6M D ^ 'As/\vork" 


chanp.e to i!ie working space 


cd SWORR 

fi.'ffO'j ipsciidopotjntials to the working space #### 


cp SPSHUDOas aa.cpi fort. 1 1 


make input file for cp-prograni; inp.mod # 
n leor Hullin'- main features of the ftii96md) # 


cal ■ inp.mod * • end B 
-1 !t)() KMHHKK) 

loot) 1 

20.0 0.9 

4.0 0.2 0.0001 
400 2 

0.0 1.0 
2 2 

0 .false. 

.F. 0.001 .F. 0.002 
.false, .false, .false. 1800 
.false. 

0.0001 0.0005 0.3 
0.001 0.001 3 

1 

end^B 


nbeg iprint timequeue 
nomore nomore_init 
delt gamma 

dclt2 gamma2 eps_chg_dlt 
delt_ion nOrder 
pff\_store mesh_accuracy 
idyn i_edyn 
i_xc t_postc 

tranc ampre tranp amprp 
tfor tdyn tsdp nstepe 
; tdipol 

; epsel epsfor epsekinc 
: force_eps niax_no_force 
; init basis 
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# 

# 

# 

# 

# 


create input file : start.ini 
(controlling the atomic structure and the input) 
and parameter.h (used for compiling of fhi96md) 
by the help of the program start.f 
make inp for start program 


# 

# 

# 

# 

# 



cat > start.inp « end_A 
1 

0 

5 

4 0 

14.2102 I 4027 0.0 1 1 1 
1 

0.33 0.33 0.25 1.0 
1 1 1 
-true. 

20 4.0 
0.1 .true. -f. 

.true, .false. 1 

5 2 1234 

873.0 1400.0 le8 1 
.t. .true. 

24 5 'Arsenic' 

1.0 74.92 0.7 3 3 
.t. .t. .f. 

0.0 0.0 4.5260 

0.0 0.0 -4.5260 

2.3684 4.7368 17.8143 

2.3684 4.7368 8.7622 

4.7368 2.3684 11.1702 

4.7368 2.3684 2.1181 

7.1051 0.0 4.5260 

7.1051 0.0 -4.5260 

9.4735 4.7368 17.8143 
9.4735 4.7368 8.7622 
11.8419 2.3684 11.1702 
11.8419 2.3684 2.1181 
0.0 7.1051 4.5260 

0.0 7.1051 -4.5260 

2.3684 11.8419 17.8143 

2.3684 11.8419 8.7622 

4.7368 9.4735 11.1702 

4.7368 9.4735 2.1181 

7.1051 7.1051 4.5260 


: number of species 
: excess electrons 
: number of empty states 
: ibrav pgind 
; celldm 

: number of k-points 
; k-point coordintes, weight 
: fold parameter 

: k-point coordinates relative or absolut? 

; Ecut [Ry], Ecuti [Ry] 

: ekt tmetal tdegen 
: tmold tband nrho 
: npos nthm nseed 
; T_ion T_init Q nfi_rescale 
: tpsmesh coordwave 
: number of atoms, zv, name 
; gauss radius, mass, danipmg,l_max,l_loc 
: t init basis 


.t. : 

tauO tford 

.t. : 

tauO tford 

.t. : 

tauO tford 

.t. : 

tauO tford 

.t. : 

tauO tford 

.t. : 

tauO tford 

.t. : 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 

.t. 

tauO tford 
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7.1051 7.1051 -4.5260 
0.4735 11.8419 17.8143 
9 4735 11.8419 8.7622 
11.8419 9.4735 11.1702 
11.8419 9.4735 2.1181 
end A 


■t. : tauO tford 

•t. : tauO tford 

•t ; tauO tford 

•t. • tauO tford 

•t. : tauO tford 




if 
P 
# 

# 

mmur 


run start program and build up parameter.h 
and inp.mi 

-in principle one could create inp.ini and 
parameter.h by hand, but start gives a 
consistent and optimized input for the 
thi96md 


# 

# 

# 

# 

# 

# 


pushd S I SOURCES} /start 
make 

if (Sstatus) then 

echo 'makes.start: compilation failed' 
exit 
endif 
popd 

cp S{ SOURCES }/start.x . 
start.x jtce start.out 




.UUUUU^JJ.44- 


a 


compile fhi96md 


# 




cp $ {SOURCES }/flii96md/parameter.h parameter.old 
set changes = 'cmp parameter.h parameter.old' 
pushd S{SOURCES}/fhi96md 

if( "Schanges" != "" ) then 
cp SWORK/parameter.h . 
make clean 


endif 

make netlib 
if (Sstatus) then 

echo 'makes.fhi96md: compilation failed' 
exit 
endif 
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popd 

m\ S {SOURCES }/flii96md.x. 


run {lii96md # 


fiii96md,x 



Appendix E 

Neural network inputs used for the optimization of arsenic 


2a (a u ) 

c (a u.) 

z 

Total Energy (a.u.) 

12 78^)3 

17 9391 

0.2044 

-148.592756 

12 7893 

17.9391 

0.2158 

-148.597763 

12 7893 

17.9391 

0.2271 

-149.685384 

12 7893 

17.9391 

0.2385 

-151 652316 

12 7893 

17.9391 

0.2498 

-151.867524 

12.7893 

18.9357 

0.2044 

-148.900903 

12 7893 

18 9357 

0.2158 

-150.334842 

12 7893 

18 9357 

0.2271 

-151335916 

12.7893 

18.9357 

0.2385 

-151.928675 

12 7893 

18.9357 

0.2498 

-152 104528 

12 7893 

19.9323 

0.2044 

-149.174143 

12.7893 

19.9323 

0 2158 

-150.476578 

12.7893 

19 9323 

0.2271 

-151.580436 

12 7893 

19.9323 

0.2385 

-152.125691 

12.7893 

19.9323 

0.2498 

-152.276857 

12 7893 

20.9289 

0.2044 

-149436043 

12.7893 

20.9289 

0.2158 

-150 881858 

12.7893 

20.9289 

0 2271 

-151.810842 

12.7893 

20.9289 

0.2385 

-152.290418 

12.7893 

20.9289 

0.2498 

-152.427439 

12 7893 

21 9255 

0.2044 

-149.405347 

12.7893 

21 9255 

0.2158 

-151.12387 

12 7893 

21.9255 

0.2271 

-151.996896 

12.7893 

21.9255 

0.2358 

-152.413284 

12 7893 

21.9255 

0.2498 

-152.532102 

13.4998 

17.9391 

0.2044 

-149.966719 

13.4998 

17.9391 

0.2158 

-150.947268 

13 4998 

17 9391 

0.2271 

-151.648953 

13 4998 

17 9391 

0.2358 

-152.034279 

13 4998 

17 9391 

0.2498 

-152.145701 
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2a (a u ) 

c (a u ) 

z 

Total Energy (a.u.) 

13 4998 

18.9357 

0.2044 

-150 190956 

13 4998 

18 9357 

0.2158 

-151 168662 

13.4998 

18 9357 

0.2271 

-151 830942 

13 4998 

18.9357 

0 2385 

-152 182843 

13 4998 

18 9357 

0 2498 

-152 297785 

13.4998 

19.9323 

0.2044 

-150.190956 

13.4998 

19.9323 

0.2158 

-151 353098 

13.4998 

19 9323 

0.2271 

-151 973077 

13 4998 

19 9323 

02385 

-152 305185 

13 4998 

19.9323 

0.2498 

-152 403676 

13.4998 

20.9289 

0.2044 

-150 549596 

13 4998 

20.9289 

02158 

-151.510053 

13.4998 

20.9289 

0.2271 

-152.093085 

13 4998 

20 9289 

0.2385 

-152.402923 

13.4998 

20.9289 

0.2498 

-152.489102 

13.4998 

21.9255 

0.2044 

-150.704569 

13.4998 

21.9255 

0.2158 

-151 655003 

13 4998 

21.9255 

0 2271 

-152.206141 

13.4998 

21 9255 

0 2385 

-152.47314 

13 4998 

21.9255 

0.2498 

-152.546772 

14.2103 

17.9391 

0.2044 

-150.896989 

14.2103 

17.9391 

0.2158 

-151.574875 

14.2103 

17.9391 

0.2271 

-152.041176 

14 2103 

17.9391 

0.2385 

-152.283256 

14.2103 

17.9391 

0.2498 

-152.3393 

14 2103 

18.9357 

0.2044 

-150.998672 

14 2103 

18.9357 

0.2158 

-151.686898 

14 2103 

18.9357 

0.2271 

-152.12534 

14.2103 

18.9357 

0.2385 

-152 338327 

14 2103 

18.9357 

0.2498 

-152 392332 



2a (a u ) 


c (a.u.) 


z 


Total Energy {a.u.) 


14.2103 
14 2103 
14 2103 

14.2103 
14 2103 

14 2103 
14 2103 
14 2103 
14 2103 

14.2103 

14.2103 

14.2103 

14.2103 
14 2103 
14 2103 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 

14.9208 


19.9323 

19 9323 

19.9323 

19.9323 

19.9323 

20.9289 

20 9289 

20 9289 

20.9289 

20.9289 

21 9255 

21.9255 

21.9255 
21 9255 

21.9255 

17.9391 

17.9391 

17.9391 

17.9391 

17.9391 

18.9357 

18.9357 

18.9357 

18.9357 

18.9357 

19.9323 

19.9323 

19.9323 

19.9323 

19.9323 


0 2044 
02158 
0 2271 
0 2385 
0.2498 

0.2044 
02158 
0.2271 
0 2385 
0 2498 

0.2044 
0.2158 
0.2271 
0.2385 
0 2498 

0.2044 

0.2158 

0.2271 

0.2385 

0.2498 

0.2044 

0.2158 

0.2271 

0.2385 

0.2498 

0.2044 

0.2158 

0.2271 

0.2385 

0.2498 


-151 190622 
- 151.837567 
- 152.214004 
- 152.382564 
- 152.423512 

- 151.300884 

- 151.92332 

- 152.263325 

- 152.417692 

- 152.463261 

- 151.381373 

- 152.036796 

- 152.309305 

- 152.426381 

- 152.458771 

- 151.507124 

- 151.96845 

- 152.267596 

- 152.410306 

- 152.430264 

- 151.618496 

- 152.062044 

- 152.318991 

- 152.423884 

- 152.425016 

- 151.705588 

- 152.125665 

- 152.339689 

- 152.407047 

- 152.403438 
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2a (a u ) 

14 9208 
14 9208 
14 9208 
14 9208 
14 9208 

14 9208 
14 9208 

14 9208 
14.9208 
14.9208 

15 6313 

15.6313 

15.6313 

15.6313 

15.6313 

15.6313 

15.6313 

15.6313 

15.6313 

15.6313 

15 6313 

15.6313 
15 6313 
15 6313 

15.6313 

15 6313 
15 6313 

15.6313 

15.6313 

15.6313 


c (au ) 

z 

Total Energy (a 

20 9289 

0.2044 

-151.774869 

20 9289 

0.2158 

-152.166093 

20 9289 

0.2271 

-152 337961 

20 9289 

0.2385 

-152.384107 

20.9289 

0.2498 

-152.387221 

21 9255 

0.2044 

-151.841269 

21 9255 

0 2158 

-152 195787 

21 9255 

0.2271 

-152 33325 

21 9255 

0.2385 

-152 364998 

21 9255 

0.2498 

-152.374562 

17 9391 

0.2044 

-151 884662 

17.9391 

02158 

-152 188452 

17.9391 

0.2271 

-152 371035 

17 9391 

0 2385 

-152 442035 

17 9391 

0.2498 

-152.443712 

18.9357 

0 2044 

-151.956486 

18.9357 

02158 

-152.23994 

18 9357 

0.2271 

-152.381797 

18 9357 

0.2385 

-152.41976 

18.9357 

0 2498 

-152 407368 

19.9323 

0,2044 

-152.008415 

19 9323 

0.2158 

-152.267179 

19.9323 

0 2271 

-152.366205 

19.9323 

0.2385 

-152.370041 

19.9323 

0.2498 

-152.351764 

20.9289 

0.2044 

-152.047868 

20 9289 

0.2158 

-152.27578 

20.9289 

0 2271 

-152.335951 

20 9289 

0 2385 

-152.313594 

20.9289 

0 2498 

-152.289382 



2a {a u ) 

c (a.u ) 

15 6313 

21 9255 

15.6313 

21 9255 

15 6313 

21 9255 

15 6313 

21 9255 

15 6313 

21.9255 


z 


0 2044 
0 2158 
0.2271 
0.2385 
0 2498 


Total Energy (a.u 

-152.078463 
-152 272436 
-152.297813 
-152.251902 
-152.227326 



Appendix F 

Input file used for the optimisation of Ga 


•’bin csh -xvf 
UUitUit?tn define home 


set FHI96MD = '’~/fhi96md" 


######??####### set directories ■ 


set PSEUDO = "S{FHI96MD}/pseudo" 
set SOURCES = ’'${FHI96MD}/src" 
set WORK = "S { FHI96MD } /Ga/work" 


mmmmum change to the ■ 


cd SWORK 


#### move pseudopotentials to the working space #### 
cp $PSEUDO/ga_aa.cpi fort.ll 


# make input file for fhi96md; inp.mod 

# (controlling main features of the fhi96md) 




# 

# 


cat > inp.mod « end_B 
-1 100 1000000 
5000 2 

30.0 0.9 

4.0 0.2 0.0001 
400 2 

0.0 1.0 
22 

0 .false. 

.F. 0.001 .F. 0.002 
.false, .false, .false. 1800 
.false. 

0.0001 0.0005 0.2 
0.001 0.001 3 
3 

end B 


: nbeg iprint timequeue 
; nomore nomorejnit 
: delt gamma 

: delt2 gainma2 eps_chg_dlt 
; delt_ion nOrder 
: pfft_store mesh_accuracy 
: idyn i_edyn 
: i_xc tjpostc 
: trane ampre tranp amprp 
; tfor tdyn tsdp nstepe 
: tdipol 

: epsel epsfor epsekinc 
: force_eps max_no_force 
; init basis 
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P create input file : inp.ini # 

U (controlling the atomic structure and the input) # 

# and parameter.h (used for compiling of fhi96md) # 

# by the help of the program start.f # 

# make inp for start program # 


cat > start, inp « end_A 


1 

; number of species 

0 

; excess electrons 

5 

; number of empty states 

8 0 

; ibrav pgind 

17.08 28.946 8.552 1 1 1 

; celldm 

1 

: number of k-points 

0.5 0.5 0.5 1.00 

: k-point coordintes, weight 

1 1 1 

: fold parameter 

.true. 

; k-point coordinates relative or absolut? 

20 5.0 

: Ecut [Ry], Ecuti [Ry] 

0.1 .true, .f 

; ekt tmetal tdegen 

.true, .false. 1 

: tmold tband nrho 

5 2 1234 

; npos nthm nseed 

873.0 1400.0 le8 1 

: T ion T_init Q nfi_rescale 

.t. .true. 

; tpsmesh coordwave 

32 3 'Gallium' 

: number of atoms, zv, name 

1.0 69.72 0.7 3 3 

: gauss radius, mass, damping, l_max,l_loc 

.t. .t. .f 

: t initjjasis 

0.0 2.2274 0.6824 

.t. : tauO tford 

0.0 -2.2274-0.6824 

.t. : tauO tford 

4.27 2.2274 3.5936 

,t. : tauO tford 

4.27 -2.2274 4.9584 

.t. : tauO tford 

4.27 9.4639 0.6824 

.t. : tauO tford 

4.27 5.0091 -0.6824 

.t. : tauO tford 

0.0 5.0091 4.9584 

.t. : tauO tford 

0.0 9.4639 3.5936 

.t. : tauO tford 

8.54 2.2274 0.6824 

.t. : tauO tford 

8.54 -2.2274-0.6824 

.t. : tauO tford 

12.81 2.2274 3.5936 

.t. ; tauO tford 

12.81 -2.2274 4.9584 

.t. : tauO tford 

12.81 9.4639 0.6824 

.t. : tauO tford 

12.81 5.0091 -0.6824 

.t. ; tauO tford 

8.54 5.0091 4.9584 

.t. : tauO tford 

8.54 9.4639 3.5936 

.t. : tauO tford 

0.0 16.7004 0.6824 

.t. : tauO tford 

0.0 12.2456-0.6824 

.t. : tauO tford 

4.27 16.7004 3.5936 

.t. : tauO tford 
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4.27 12.2456 4.9584 

.t. 

: tauO tford 

4.27 23.9369 0.6824 

.t. 

: tauO tford 

4.27 19.4821 -0.6824 

.t. 

: tauO tford 

0.0 19.4821 4.9584 

.t. 

: tauO tford 

0.0 23.9369 3.5936 

•t. 

: tauO tford 

8.54 16.7004 0.6824 

.t. 

: tauO tford 

8.54 12.2456-0.6824 

.t. 

: tauO tford 

12.81 16 7004 3.5936 

.t. 

: tauO tford 

12.81 12.2456 4.9584 

.t. 

; tauO tford 

12.81 23.9369 0.6824 

■t. 

: tauO tford 

12 81 19 4821 -0.6824 

.t. 

: tauO tford 

8.54 19.4821 4.9584 

.t. 

: tauO tford 

8.54 23.9369 3.5936 
end A 

.t. 

• tauO tford 




# run start program and build up parameter.h # 

U and inp.mi # 

tf -in principle one could create inp.ini and # 

# parameter.h by hand, but start gives a # 

# consistent and optimized input for # 

# fhi96md # 


pushd ${ SOURCES }/start 
make 

if (Sstatus) then 

echo 'makes. start: compilation failed' 
exit 
end if 
popd 

cp ${ SOURCES }/start.x . 
start.x |tec start. out 



cp ${SOURCES}/fhi96md/parameter.h parameter.old 
set changes = emp parameter.h parameter.old 
pushd ${SOURCES}/fhi96md 

if( "Schanges" != "" ) then 
cp SWORK/parameter.h . 
make clean 
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endif 


make netlib 
if (Sstatus) then 

echo 'makes. fl'ii96md. compilation failed' 
exit 
endif 
popd 

mv S{ SOURCES }/fhi96md.x . 


# run fhi96md # 

###################################################### 


flii96md.x 


103 



Appendix G 

Neural network inputs used for the optimization of gallium 


2a (a u.) 

2b (a u ) 

c (a u.) 

y 

z 

Total energy (a u ) 

15 372 

15 372 

26 0507 

26.0507 

7 6973 

7 6973 

0 1385 

0.1385 

0.0718 

0.0878 

-72 946632 

-73.00116 

15 372 

15 372 

26 0507 

26 0507 

7 6973 

7.6973 

0 1693 

0 1693 

0.0718 

0.0878 

-73 17925 

-73 116578 

15 372 

15.372 

26 0507 

26.0507 

8.5525 

8 5525 

0.1385 

0.1385 

0.0718 

0.0878 

-72 717299 

-72 802232 

15.372 

15 372 

26.0507 

26.0507 

8 5525 

8 5525 

0.1693 

0 1693 

0 0718 

0.0878 

-73 042343 

-73.01291 

15.372 

15.372 

26 0507 

26.0507 

9.4078 

9 4078 

0 1385 

0.1385 

0.0718 

0.0718 

-72 325745 

-72 4347 

15 372 

15.372 

26 0507 

26.0507 

9 4078 

9 4078 

0 1693 

0 1693 

0 0878 

0.0878 

-72 649952 

-72 630669 

15.372 

15.372 

28.9452 

28.9452 

7 6973 

7.6973 

0.1385 

0.1385 

0 0718 

0 0878 

-73 282178 

-73.274833 

15.372 

15.372 

28.9452 

28.9452 

7 6973 

7.6973 

0.1693 

0.1693 

0.0718 

0.0878 

-73 301977 

-73 25467 

15.372 

15.372 

28.9452 

28.9452 

8 5525 

8 5525 

0 1385 

0.1385 

0.0718 

0 0878 

-73 104797 

-73 108821 

15.372 

15.372 

28.9452 

28.9452 

8 5525 

8.5525 

0.1693 

0 1693 

0 0718 

0 0878 

-73 149786 

-73.119617 

15.372 

15.372 

28.9452 

28.9452 

9 4078 

9 4078 

01385 

0 1385 

0 0718 

0 0878 

-72 650385 

-72 686652 

15 372 

15.372 

28.9452 

28.9452 

9.4078 

9 4078 

0.1693 

0.1693 

0 0718 

0 0878 

-72 764318 

-72 748086 



2d id 12 1 

’b la u ) 

(a u ) 

1 3 

1 3"; 

.N3‘r 

'! S3>r 

7 6973 
7 6973 

I *3 

15 3"? 

'] Siir 

M 

^’.6973 
7 6973 

1 5 3 

15 3 T 2 

'! S5‘)'' 

.'1 S3‘)T 

S 5525 
8 5525 

is 

1 5 3 

?1 S?‘)- 

M S.?')” 

8 5525 
8 5525 

1^ V*: 

I ^ ' ‘3 

' 1 S.U)-’ 

M 

9.4078 

9.4078 

l> ^ '3 

1 > ' "3 

'' ! S 

M ,s;‘r 

9.4078 

9 4078 

!'■ 118 

n ii8 

SKir 

.■’<>()S(I7 

7 6973 

7 6973 

!7Jt8 

!7(KS 

26 ()5cr 

26 0507 

7 6973 
7,6973 

17 08 

17 08 

26,()7()7 

26 0507 

K.5525 

8,5525 

17.(18 

I 3.08 

2<. OSO" 

26 070? 

8.5525 

K 5525 

17 08 

17 08 

26 U>()7 

26,0507 

9.4078 

9 4078 

17 08 

17 08 

26.0507 

26 0507 

9.4078 

9.4078 

17 08 

17.08 

2.S ‘>452 

2K 9452 

7.6973 

7.6973 


y 

2 

Total Energy (a. 

0 1385 

0 1385 

0.0718 

0.0878 

-73 210365 
-73.20457 

0,1693 

0 1693 

00718 

0 0878 

-73.280907 
-73 248118 

0.1385 

0.1385 

0 0718 
0.0878 

-73.013472 

-73.029739 

01693 

0.1693 

0 0718 

0 0878 

-73 126805 
-73.112496 

0.1385 

0 1385 

0 0718 

0 0878 

-72.646568 
-72 681739 

0.1693 

0.1693 

0.0718 

0 0878 

-72 759612 
-72 762634 

0.1385 

0.1385 

0 0718 

0 0878 

-72 87244 
-72 95765 

0 1693 
0.1693 

0.0718 

0 0878 

-73 174164 
-73 137081 

0.1385 

0.1385 

0.0718 

0.0878 

-72.670453 

-72.775988 

0.1693 

0.1693 

0.0718 

0.0878 

-72.067294 

-73.037114 

0.1385 

0.1385 

0.0718 

0.0878 

-72.312578 

-72.431531 

0.1693 

0 1693 

0.0718 

0 0878 

-72.706837 
-72 674234 

0 1385 

0.1385 

0.0718 

0 0878 

-73.152382 

-73.163748 



2a (a.u.) 

2b (a.u.) 

c (au) 

y 

z 

Total Energy (a u ) 

17.08 

17 08 

28.9452 

28.9452 

7.6973 

7.6973 

0.1693 

0.1693 

0.0718 

0 0878 

-73 198104 

-73 153634 

17 08 

17 08 

28 9452 

28 9452 

8.5525 

8.5525 

0.1385 

0.1385 

0 0718 

0 0878 

-72 927961 

-72 95572 

17.08 

17.08 

28.9452 

28 9452 

8.5525 

8 5525 

0.1693 

0.1693 

0.0718 

0.0878 

-73 077214 

-73 038365 

17.08 

17 08 

28.9452 

28.9452 

9.4078 

9.4078 

0.1385 

0.1385 

0.0718 

0.0878 

-72 580664 

-72 613028 

17.08 

17.08 

28.9452 

28.9452 

9.4078 

9.4078 

0.1693 

0 1693 

0.0718 

0.0878 

-72.721315 

-72 685313 

17.08 

17.08 

31.8397 

31.8397 

7.6973 

7 6973 

0 1385 

0.1385 

0.0718 

0.0878 

-73 023889 

-73.027779 

17.08 

17 08 

31.8397 

31.8397 

7.6973 

7.6973 

0 1693 

01693 

0.0718 

0.0878 

-73.095547 

-73 065969 

17 08 

17 08 

31.8397 

31.8397 

8.5525 

8.5525 

0 1385 

0.1385 

0.0718 

0.0878 

-72 854404 

-72 869299 

17 08 

17.08 

31.8397 

31.8397 

8.5525 

8.5525 

01693 

0.1693 

0.0718 

0.0878 

-72 9741 

-73.95221 

17 08 

17.08 

31.8397 

31.8397 

9.4078 

9.4078 

0.1385 

01385 

0.0718 

0 0878 

-72 521119 

-72 543919 

17 08 

17 08 

31.8397 

31.8397 

9.4078 

9.4078 

01693 

01693 

0.0718 

0.0878 

-72 637968 

-72.623488 

18.788 

18.788 

26.0507 

26.0507 

7.6973 

7.6973 

01385 

0 1385 

0.0718 

0.0878 

-72 653423 

-72.754084 

18.788 

18.788 

26.0507 

26.0507 

7.6973 

7.6973 

0 1693 

01693 

0 0718 

0 0878 

-72 99646 

-72 943921 



2a (a.u.) 

2b (a.u.) 

c (a.u.) 

y 

z 

Total Energy (a u ) 

18 788 

18.788 

26.0507 

26 0507 

8.5525 

8.5525 

0.1385 

0.1385 

0 0718 

0 0878 

-72.478426 

-72 595728 

18 788 

18 788 

26.0507 

26 0507 

8 5525 

8.5525 

0 1693 

0.1693 

0 0718 

0 0878 

-72.916775 

-72 88746 

18 788 

18 788 

26 0507 

26.0507 

9.4078 

9 4078 

0 1385 

0.1385 

0.0718 

0.0878 

-72.152547 

-72.279387 

18 788 

18 788 

26 0507 

26 0507 

9.4078 

9.4078 

0 1693 

0 1693 

0 0718 

0 0878 

-72 581597 

-72 548741 

18 788 

18 788 

28 9452 

28.9452 

7.6973 

7.6973 

0.1385 

0.1385 

0 0718 

0 0878 

-72 918081 

-72.937145 

18 788 

18 788 

28.9452 

28 9452 

7.6973 

7.6973 

0.1693 

0.1693 

0 0718 

0.0878 

-72.969414 

-72.926379 

18.788 

18 788 

28 9452 

28 9452 

8 5525 

8.5525 

0.1385 

0 1385 

0.0718 

0 0878 

-72 796325 

-72.810759 

18 788 

18 788 

28 9452 

28.9452 

8 5525 

8.5525 

0.1693 

0 1693 

0.0718 

0.0878 

-72.875041 

-72.832872 

18 788 

18.788 

28 9452 

28.9452 

9.4078 

9 4078 

0 1385 

0 1385 

0.0718 

0.0878 

-72 410677 

-72 441576 

18 788 

18 788 

28 9452 

28.9452 

9.4078 

9 4078 

0.1693 

0 1693 

0.0718 

0.0878 

-72.548506 

-72 504717 

18 788 

18.788 

31 8397 

31 8397 

7.6973 

7.6973 

0 1385 

0 1385 

0.0718 

0 0878 

-72 760679 

-72 758218 

18 788 

18 788 

31.8397 

31.8397 

7 6973 

7.6973 

0 1693 

0.1693 

0.0718 

0.0878 

-72 802653 

-72 767385 

18.788 

18.788 

31.8397 

31 8397 

8 5525 

8.5525 

0 1385 

0.1385 

0.0718 

0.0878 

-72.623396 

-72 625656 



-a(au) 

2b(a u ) 

c (a u.) 

y 

z 

Total Energy (a.u ) 

18.788 

18.788 

31 8397 

31 8397 

8 5525 

8 5525 

0.1693 

0.1693 

0 0718 

0 0878 

-72 710581 

-72.677826 

18 788 
18.788 

31 8397 

31 8397 

9.4078 

9.4078 

0 1385 

01385 

0.0718 

0.0878 

-72.328431 

-72.331826 

IS 788 

18 788 

31 8397 

31 8397 

9.4078 

9.4078 

0.1693 

0.1693 

0.0718 

0.0878 

-72.405067 

-72.374261 
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